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A  SIGNAL  DETECTION  FROM  NOISE 
UTILIZING  ZERO-CROSSING  INFORMATION 


By 
Woon  Cheon  Yeo 
December,  1975 

Chairman:  Dr.  Jack  R.  Smith 

Major  Department:  Electrical  Engineering 

A  method  of  detecting  signals  from  noise  utilizing  the 
zero-crossing  information  is  analyzed. 

When  the  noise  is  stationary  white  Gaussian,  a  corre- 
lation detector  is  considered  to  be  optimal.  But,  if  the 
process  of  signal  and  noise  is  complex,  then  a  nonlinear 
detector  may  give  a  superior  result  in  overall  performance. 

In  a  situation  where  the  signals  and  noise  have 
nonstationary  statistics,  a  detector  which  measures  the 
zero-crossing  intervals  of  the  process  has  been  found  to  be 
very  efficient  in  detecting  signals  from  noise.  The  operating 
characteristic  of  the  zero-crossing  measuring  detector  is 
analyzed  in  a  model  with  signals  of  sine  wave  form  imbedded 
in  a  stationary  bandlimited  white  Gaussian  noise.  Then,  to 


observe  the  operating  characteristic  of  the  detector  in  a 
complex  situation,  the  detector  is  analyzed  under  deviated 
conditions  of  signals  and  noise  from  the  assumptions.  The 
performance  of  a  correlation  detector  is  also  analyzed  under 
identical  conditions,  as  a  comparative  reference  to  the 
performance  of  the  zero-crossing  measuring  detector. 

The  operating  characteristic  of  the  zero-crossing 
measuring  detector  is'  shown  to  he  effective  in  detecting 
phasic  events  from  a  sleep  electroencephalogram,  where  the 
process  of  the  signals  and  noise  is  very  complex. 


vi 


CHAPTER  I 
INTRODUCTION 


An  optimum  detector  of  a  signal  from  noise  may  well  be 
described  as  "a  detector  which  best  satisfies  given  criteria 
under  a  given  set  of  assumptions"   (Whalen). 

The  first  problem  in  designing  an  optimum  detector  is, 
then,  making  proper  assumptions  on  parameters  that  completely 
describe  the  process  of  signal  and  noise,  and  making  proper 
selections  of  parameters  that  are  relevant  to  discriminating 
signals  from  noise.  If  either  the  real  process  is  different 
from  that  described  by  the  assumptions  on  parameters,  or  any 
of  the  relevant  parameters  are  neglected,  then  the  detector 
may  be  only  an  approximation  to  the  optimum. 

In  the  case  when  the  process  is  complex,  we  may  not  be 
able  to  make  assumptions  on  parameters,  or  even  if  we  could 
define  the  process  by  the  parameters,  an  optimum  detector 
based  on  all  the  relevant  parameters  becomes  too  complicated 
to  implement.  Thus  we  need  to  select  among  the  definable 
parameters  the  most  efficient  featiires  from  the  process  as 
the  discriminating  parameters  so  that  a  detector  based  on 
these  parameters  should  give  an  acceptable  performance  in 
the  complex  situation.  The  efficiency  of  the  parameters  is 
judged  by  the  given  criteria  under  the  assumptions  on  the 
parameters.  The  performance  of  the  detector  based  on  these 


efficient  parameters,  however,  may  not  be  optimum  if  the 
environment  of  signal  and  noise  changes  so  as  to  include 
other  efficient  definable  parameters. 

In  a  stationary  Gaussian  noise,  a  correlation  detector 
is  considered  to  be  optimal  (Wainstein) .  If  the  noise  process 
is  complex  so  that  it  could  not  be  simply  definable  as  sta- 
tionary Gaussian,  then  a  superior  detector  may  require  a 
nonlinear  characteristic,  and  the  implementation  of  the 
detector  may  be  complicated  and  costly.  Especially  if  the 
process  is  not  stationary,  then  the  adaptation  of  the  detec- 
tor becomes  a  serious  problem,  and  it  may  leave  the  detector 
complicated  to  operate.  Thus,  in  general,  the  design  of  a 
detector  should  include  considerations  such  as  performance, 
cost  of  implementation,  simplicity  of  implementation,  and 
simplicity  of  operation.  The  priority  of  the  considerations 
may  vary  depending  on  the  circumstances  of  the  particular 
application  of  the  detector. 

In  designing  a  system  to  detect  phasic  events  from 
sleep  electroencephalograms (EEGs ) ,  one  encounters  the 
following  complexities  in  the  process  of  the  EEG  activities: 
The  phasic  events,  which  are  designated  by  arrows  in  Figure 
23 »  are  distinguished  from,  the  irregular  background  EEG  by 
their  resemblance  to  sinusoidal  waves.  Sigma  spindles,  for 
example,  last  from  6  to  20  cycles  with  a  frequency  between 
12  and  15  Hz  depending  on  the  individual.  The  amplitude  of 
spindles  and  background  EEG  activities,  which  is  considered 
as  noise  when  we  desire  to  detect  sigma  spindles,  also 


varies  widely  depending  on  the  individual  and  the  electrode 
resistance.  The  a  priori  probability  densities  of  the  fre- 
quency and  amplitude  distribution  are  not  usually  known. 
Even  though  the  strength  of  the  signal  (sigma  spindles)  and 
noise  varies,  the  signal-to-noise  ratio  remains  relatively 
constant;  that  is,  a  subject  with  high  background  EEG  acti- 
vity also  has  high  amplitude  in  the  spindle  activity.  The 
strength  and  the  power  spectra  of  the  background  EEG  undergo 
a  continuous  change  as  a  sleep  progresses. 

In  practical  applications  of  the  phasic  event  detector, 
it  is  desirable  that  the  detector  be  able  to  perform  equally 
well  an  any  subject  without  a  priori  information  on  the 
frequency  and  the  strength  of  the  EEG  activity.  To  accomplish 
this  a  detector  could  be  designed  which  adapts  to  the  varying 
frequency  and  intensity  of  the  EEG  activity,  but  it  might  be 
too  costly  to  implement. 

Smith  et  al, {k3,kk)   designed  a  phasic  event  detector 
which  employs  the  zero-crossing  interval  of  the  EEG  as  the 
discriminating  parameter,  and  they  obtained  rather  superior 
performance  from  the  detector.  The  detection  rate  was  high 
enough  to  meet  the  requirement,  while  the  false  alarm  rate 
was  extremely  low  for  any  kind  of  noise  situation.  The 
detector  could  be  applicable  to  a  wide  variety  of  subjects 
without  any  adjustment  of  thresholds  to  accommodate  each 
subject.  It  was  simple  to  implement,  and  its  overall  perfor- 
mance was  superior  to  a  correlation  detector  based  on  weak 
assumptions  on  the  signal  and  noise.  And  thus  a  mathematical 


analysis  of  the  zero-crossing  measuring  detector  v/as 
suggested  and  carried  out  in  this  study  to  provide  a  better 
understanding  on  the  efficiency  of  the  feature  of  zero- 
crossing  interval  as  a  discriminating  parameter  of  signals 
from  noise. 

Since  both  the  signal(phasic  events)  and  the  noise(back- 
ground  EEC-  activity)  processes  in  sleep  EEG  are  too  compli- 
cated to  be  explicitly  defined,  an  analysis  of  any  phasic 
event  detector  from  sleep  EEG  is  necessarily  based  on  improper 
assumptions.  Also  any  attempt  to  make  the  model  close  to  the 
reality  would  leave  the  analysis  untenable.  Thus  the  analysis 
of  the  detector  is  carried  on  with  an  idealized  model  of 
signal  and  noise.  The  signal  is  modeled  by  sinusoidal  waves 
and  the  noise  is  assumed  to  be  bandlimi ted  stationary  Gaussian. 
Then,  to  observe  the  performance  of  the  detector  in  a  nonideal 
situation,  the  operating  characteristic  of  the  detector  is 
analyzed  under  conditions  of  signal  and  noise  different  from 
the  assumptions.  In  each  case  the  performance  of  the  correl- 
ation detector  is  also  analyzed  as  a  comparative  reference  to 
the  zero-crossing  measuring  detector. 

To  compute  the  detection  and  false  alarm  probabilities  of 
the  zero-crossing  measuring  detector,  we  need  the  zero-cross- 
ing interval  density  functions  of  the  processes  with  noise 
alone  and  signal  plus  noise.  The  distribution  functions  for 
the  successive  zero-crossing  intervals  of  the  bandlimited 
Gaussian  noise  have  been  developed  by  Rice(36,37),  McFadden 
(26,2?)  and  Longuet-Higgins(23,24).  Experimental  works  on 


the  zero-crossing  interval  distribution  of  Gaussian  noise 
have  been  presented  "by  Rainal(32)  and  Blotekjaer{6) . 
Rice's  zero-crossing  inteirval  density  function  agrees  well 
with  experimental  data  for  small  values  of  zero-crossing 
interval.  The  zero-crossing  interval  density  function  derived 
by  McFadden  gives  good  approximations  to  Gaussian  processes 
having  arbitrary  power  spectra.  For  the  particular  Gaussian 
process  which  has  an  ideal  low  pass  spectrum,  the  approxi- 
mation by  Longuet-Higgins  fits  closely  to  the  experimental 
data  given  by  Blotekjaer.  Longuet-Higgins  derived  a  proba- 
bility density  function  of  the  spacing  between  the  ith  zero 
and  the  (i+m+l)th  zero  of  a  stationary  random  function.  The 
density  function  is  expressed  as  a  rapidly  convergent  series, 
and,  when  the  process  of  the  random  function  is  Gaussian,  the 
first  two  terms  of  the  series  may  be  expressed  in  terms  of 
known  functions.  This  Longuet-Higgins  approximation  is  used 
here  to  compute  the  zero-crossing  interval  distribution  of  a 
bandlimited  Gaussian  process. 

In  the  signal-plus-noise  case,  the  zero-crossing  inter- 
val density  function  is  approximated  by  the  first  term  of 
the  Longuet-Higgins^  series.  This  approximation  is  also 
derived  by  Rice (38)  for  Gaussian  processes.  Cobb (10)  gave  a 
treatment  on  this  approximation  for  a  process  of  sine  wave 
plus  bandlimited  Gaussian  noise.  He  reports  in  the  same  paper 
about  the  zero-crossing  interval  measioring  principle  in 
distinguishing  two  signals  separated  in  frequency.  Yet  his 
paper  does  not  provide  details  of  the  detector  and  any 


comparative  result  to  other  kinds  of  detecting  methods.  The 
coraputation  of  the  zero-crossing  interval  density  for  sine 
wave  plus  noise  in  this  study  is  based  on  the  treatment  of 
Cobb. 

Derivations  of  the  zero-crossing  interval  density  func- 
tion for  Gaussian  noise  and  Gaussian  noise  plus  signal  are 
presented  in  Chapter  II,  An  expression  of  the  first  two  terms 
of  the  Longuet-Higgins*  series  in  terms  of  the  autocorrelation 
function  of  the  Gaussian  noise  is  derived.  The  computed 
density  function  is  plotted  along  with  the  experimental  data 
as  a  comparison.  The  derivation  and  the  computed  data  of  the 
zero-crossing  interval  density  for  signal  plus  noise  are  also 
presented  in  Chapter  II. 

In  Chapter  III,  the  schematic  of  the  zero-crossing 
measuring  detector  is  described  and  the  method  of  maximizing 
the  detection  probability  for  given  false  alarm  rate  is 
discussed.  The  performance  of  the  detector  is  compared  to 
that  of  the  conventional  correlation  detector  under  the  iden- 
tical situation  of  signal  and  noise . 

The  advantages  and  disadvantages  of  the  zero-crossing 
measuring  detector  compared  to  the  correlation  detector  are 
further  discussed  in  Chapter  IV,  with  practical  considerations 
in  detecting  phasic  events  from  sleep  EEGs. 

All  the  programs  to  compute  various  functions  were 
written  in  FORTRAN  and  assembly  language,  and  they  were 
processed  through  a  minicomputer  PDP-8/e.  Without  an  extended 
ar5,thmetic  unit,  the  speed  of  the  PDP-8/e  was  very  slow  to 
compute  some  numerical  integrations.  Thus »  to  save  the 
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computing  time,  we  employed  data  interpolations  in  the 
computation  of  the  zero-crossing  interval  density  functions 
and  evaluated  the  statistical  functions  from  their  series 
expansions. 


CHAPTER  II 
ZERO -CROSSING  INTERVAL  DENSITY  FUNCTIONS 


Derivation  of  Zero-crossing:  Interval  Density  Function 

For  simplicity,  we  will  state  that  a  random  function 
f(t)  has  a  zero  or  a  crossing  at  t=t.  if  f(t. )=0.  Similarly 
we  use  up-crossing  or  down-crossing  to  indicate  a  crossing 
with  f'(t^)>0  or  f'(t^)<0  respectively. 

Consider  the  instance  when  a  random  function  f (t)  has 
up-crossings  in  the  infinitesimal  intervals  (t.,t^+dt.) 
(i=l,.,,,n),  and  let  us  denote  the  joint  probability  of  this 
happening  by  W(+,+,  ...,+)dt^dt2.  .  .dt^^.  To  denote  the  joint 
probability  that  the  random  function. has  down-crossings  in 
some  intervals,  we  change  the  plus  sign   to  minus  at  the 
corresponding  positions  in  W.  Thus  W(+, -,-}-)dt.  dt^dt^,  for 
instance,  would  represent  the  joint  probability  that  a  random 
function  f(t)  has  an  up-crossing  in  the  interval  (t. ,t^+dt. ), 
dovm-crossing  in  (t^, t^+dt^) ,  and  an  up -crossing  in 
(t^,t^+dt^). 

V/e  next  introduce  pCmjr)  to  denote  the  probability 
density  of  the  interval  between  the  ith  and  the  (i+m+l)th 

zeros  of  f(t).  This  probability  density  is  denoted  by  P  (r) 

m   ' 

by  other  authors,  but,  since  in  this  paper  a  probability  is 
represented  by  an  uppercase  P  and  a  probab5.1ity  density  is 


9 

v/ritten  in  a  lowercase  p,  the  symbol  p(in;r)  is  used  here. 
Then  p(0;r)d'r  is  the  probability  that  the  interval  between 
two  successive  zeros  has  length  t»  and  p(0;r)drW(+)dt^ 
expresses  the  probability  that  f (t)  has  an  up-crossing  at 
t=t.  and  a  down-crossing  after  an  interval  -r  without  having 
any  other  crossings  in  between  the  two  crossings.  Likewise 
p(2n;r)drW(+)dt.  represents  the  probability  of  f(t)  having 
an  up-crossing  at  t=t:  and  a  down-crossing  after  an  interval 
-Ti  containing  2n  crossings  anywhere  in  between  the  two 
crossings. 

If  f (t)  has  an  up-crossing  at  t=t^  and  a  down-crossing 
at  t^tp,  then  f(t)  may  contain  2n  (n=0,l,2, . . . )  crossings 
between  t=t.  and  t=t2.  W(+, -)dt. dt^i  the  probability  of  f(t) 
having  an  up-crossing  at  t=t.  and  a  down-crossing  at  t=t2, 
is  then  the  svaa.   of  the  probabilities  of  all  possible  cases 
when  f(t)  has  an  up-crossing  at  t=t^  and  a  down-crossing 
after  an  interval  r=t2-t.,  including  2n  (n=0,l,2, . . . )  cross- 
ings  within  the  interval  r .  That  is 

W(+,-)dt^dt2=[p(0;r)+p(2;r)+p(^;r)+...]dtW(+)dt^,  (1.1) 

or,  since  dT=dtp, 

W(  +  ,-.)/W(  +  )=p(0;T)+p(2;r)+p(4;r)  +  ...,  (1.2) 

(McFadden  1958). 

In  a  similar  fashion  we  can  derive  an  expression 
related  to  W(+,-,-).  Let  us  consider  the  situation  when  f(t) 
has  an  up-crossing  at  t=t^  and  a  down-crossing  at  t=toi 
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containing  a  down-crossing  at  t=t2  anyvihere   in  between  t-  and 
to.  Since  the  crossings  at  t=t-  and  t=t~  are  an  up-crossing 
and  a  down-crossing  respectively,  f(t)  caji  only  have  an  even 
number  of  crossings  in  the  interval  (t^.t-,),  and  one  of  the 
down- crossings  among  these  crossings  should  pass  zero  at  t=t2. 
Suppose  we  have  2n  crossings  in  (t^,t-)  where  t--t^=T.  Then 
the  down-crossing  at  t=t2  could  be  one  of  those  n  down- 
crossings,  and,  therefore,  there  are  n  different  possibilities 
of  f(t)  having  a  down-crossing  at  t=t2.  Since  the  down-cross- 
ing at  t=t2  could  be  anywhere  in  (tj^,to).  if  we  integrate  the 
probability  of  each  possibility  with  respect  to  tg  from  t^  to 
to  we  get  p(2n;r)dTW(+)dt^,  and  we  obtain  for  the  n  possibili- 
ties a  total  probability  of  np(2n;r)d'''W(  +  )dt^.  And  the  sum  of 
the  probabilities  for  all  possible  values  of  n, 

[p(2;r)+2p(4;r)+...+np(2n;r)+...]drW(+)dt^,        (1.3) 

is  then  the  probability  of  f(t)  having  an  up-crossing  at  t=t^ 
and  a  down-crossing  at  t=t.+T  containing  at  least  one  down- 
crossing  at  t=t2  anywhere  in  between  t^  and  t^. 

Now  if  we  integrate  W(+,-, -)dt. dt2dto  with  respect  to  tg 
over  the  interval  (t^,to)i  we  obtain  the  probability  of  f(t) 
having  an  up- crossing  at  t=t^  and  a  down- crossing  at  t=to 
containing  at  least  one  down-crossing  at  t=t2  anywhere  in 
between  t^  and  t-.  This  probability  should  be  equal  to  that 
expressed  by  Equation  (1.3).  Therefore  we  obtain  (Longuet- 
Higgins ) , 

S        dt2[W(+,-,-)A(+)]  = 
V^2<b      p(2;T)+2p(4;r)  +  ...+np(2n;T)  +  ...        (1.^) 
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By  rearranging  Equa.tions  (1.2)  and  (l.^-), 

p(Ojr)=w(+,-)/w(+)-[p(2;r)+p(45T)+. ..],        (1.5) 

p(2jr)-  S         cltpCW(+,-,-)/W(+)]  - 

[2p(i]-;r)+3p(6;r)+...].      (1.6) 

Substituting  Equation  (1.6)  into  Equation  (1.5). 

p(0;r)=W(+,-)/W(+)  -    S       dtJLVa+»-,-)AH+)l   + 

W^3 
p(^;-r)+2p(6;r)+3p(8;-r)+....  (1.7) 

If  we  neglect  the  terms  [p(^;r)+2p(6;  t)+3p(8;t)+. . .]  in  the 
above  equation,  we  get  an  approximate  expression  for  the 
zero-crossing  interval  density  (Longuet-Higgins) : 

p(0;r)=W(  +  ,-)/W(+)-    S        dt^Cw(  +  ,-,-)/W(-J-)].    (1.8) 
t,<t2<t3 

By  neglecting  the  remaining  terms  in  Equation  (1.?)  we 
are  ignoring  the  probabilities  of  f  (t)  having  ^■   or  more 
zeros  within  the  interval  of  length  r. 
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Evaluation  of  W  for  Gaussian  Process 


Let  us  consider  the  probability  W(+,+, . . . ,+)dt. dtp. . .dt 

that  f(t)  has  an  up-crossing  with  an  arbitrary  positive  slope 
f ' (t)  in  each  of  the  small  intervals  (t. ,t^+dt^)  (i=l,...,n). 


For  convenience  write 


f(t.)=f.. 


and  let  pCfn  ,  .  . .  tfj^»g^  .  .  • .  fg„)  denote  the  joint  probability 
density  of  the  f.  and  g.  (i=l, . . . ,n) .  Thus  p(f^,...,f^, 
g^ , . . .  ,gj^)df^  .  . .  dfj^dg^  .  . .  dg^  is  the  probability  that  the  f  ^ 
and  g.  lie  in  the  intervals  (f.,f.+df.),  (g-,g^+dg^)  respec- 
tively for  each  i. 

If  f(t)  has  a  zero  in  (t.,t.+dt.)  with  a  gradient  g., 
then  f (t. )  must  lie  in  a  small  range  of  value  of  extent 
jg.|  dt. .  Especially  if  f (t)  has  an  up-crossing  in  the 
interval  (t. ,t.+dt.),  then  the  range  of  the  gradient  is 
0<g.<oo,  and  the  range  of  f.  is  -g.dt.<f.<0.  Thus  the  integral 
of  p(f^,  ...  ,f^,g^,  . .  .  ,g^)  over  the  ranges  (0<g^<°o), 
(-g^dt^<f^<0)  (i=l,...,n)  is  the  probability  of  f(t)  having 
up-crossings  at  the  intervals  (t., t-+dt.)  (i=l,...,n).  That 
is, 

W(+ +)dt^...dt^ 

OO  00  0  0 

=  J'dg.  ../dg    S     df.  ...   J*  df  p(f f  s  2) 

0  1   0  "-g^dt,  ^   -p;  dt   ^   1 in'Si""gn^- 

11       ^n  n 

(2.1) 
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we 


Since  the  range  of  f.  (-g.dt.,0)  is  small  for  each  i, 
may  consider  p(f. , . . . ,f  ,g^ g  )  to  be  constant  with 


respect  to  f.  within  each  interval  (t.,t.+dt.).  Then  the 
integrals  in  Equation  (2.1)  could  be  approximated  as 


W(+,... ,+)dt^...dt^ 

OO  00  0  0 

=/dg, .../dg  p(f f  ,g  ,...,g  )    X  df  ...    S   df  , 

0   -    0  ^,   ^      "^      ^  -g^dt^  1    -g^dt^  ^  . 

which  becomes 


W(  + +)dt^...dt^ 


and  thus , 

W(+, ...,+) 


=.rdg^...;dg.^g^...g^p(0,.  ..,0,g^,...,g^).        (2.2) 

Nov/  we  try  to  find  the  expression  of  p(f.  ,...,f  , 
gi>'..»g  )  in  terms  of  the  autocorrelation  function  of  the 
Gaussian  process  f(t).  The  covariance  matrix  of  the  2n 
variables  f ^ , . . . , f ^, g^ ^n  ^^ 


M= 


Ef;f^,....Ef^f^,Ef^g.^,...,Ef;g^ 
^^ih Eg^f^.Eg^g^....,Eg^g^_ 


(2.3) 


,  Eg^/^ .... .Eg^f^.Eg^g^ , . . . .Eg^g^  J    , 
v/here  Ef.g.  represents  Erf.g.],    an  expectation  of  f-g.. 

-1-      J  -*-      d  ■'■0 


1^      . 

Let  R. .  denote  the  autocorrelation  function  of  f (t) 
such  as 

where -r=t.-t..  Then, 

E[fiSp=E[f(t.)^f(t)|^^^  ] 

=E[f(t.)^f(t+r)jt=t.] 
=E[f(t.)^(t.+TOi  . 
=g-|E[f(t^)f(t.+T)] 

=6-|R(T)=R..'('r),  (2.5) 

where  the  prime  denotes  differentiation  with  respect  to  r, 
and  6  denotes  partial  differentiation. 


ECei^o>Eyf(*)|t=t,*(*o)^ 


=ECg^f(t-r)|^^^  f(t.)] 

J 

=-,-|E[f(t.-r)f(t,)] 


=-j|R(-T)=-5|R(T)=-Rj^..(r),  (2.6) 

1  .1 
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=E[5|f(t)j,^,,j|f(t+r)|^,.,_] 

»       1 

=J-^C-5|E[f(ti)f(t^+r)]] 

=.^^:R(r)=-R,...(r). 

Thus  the  covariance  matrix  in  terms  of  R(t)  is 


(2.7) 


^-   '  ^11' ^In'  ^ll'"-"  ^In' 


^11  '••'  ^^In  •  ^11  "••'  "in 


-^nl ' '  •  "  -^nn' '  "\l" "^nn" 


(2.8) 


By  the  Gaussian  hypothesis  we  have 
P(fi.---.fn'Si"-"Sn) 

=  (2'rr )  ""D"*exp  ( -iX^M'^X ) 


2n 


=  (2Tr)'"^^''^exp(-A  L     q.  .f .  f  . ) , 


(2.9) 


where  X  and  X-^  are  colvimn  and  row  vectors  of  2n  variables 

(f  3^ ^n'  ^1 '  •  •  • ' ^n  ^  respectively ,  D=  (M  | ,  f^+i^^i '  ^^^  °-i  j 

is  the  element  of  the  ith  row  and  .jth  column  of  the  matrix  . 
Q=M"  .  We  will  denote  a  matrix  with  elements  q.  .  by  (q^^^)' 
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Substituting  Equation  (2.9)  into  Equation  (2.2)  we  get 


W(+,...,+) 


=  (2Tr)-"D--j;dg^...j;dg^g  g2...g^exp(-i  ^  q^^.^^^.g.g.). 

(2.10) 

The  summation  in  Equation  (2.10)  involves  only  the 
last  n  rows "and  columns  of  (q^^).  We  let  (h. .)  denote  the 
inverse  matrix  of  the  matrix  composed  of  the  last  n  rows 
and  columns  of  (q^^)  such  as 

%+l ,  n+1 ^n+l ,  2n  Y  "■'■ 

(2.11) 

^2n,n+l""'V'^2n,2n     , 

and  modify  Equation  (2.10)  using  this  matrix. 


(h,3)= 


W(+ +) 


,-tn 


=  (2,r)-^"D-^|(h.  .)|  ^  Xdgj^. .  ./dg^g^gj. .  .g„[(2„)- 
|(h..)|-*exp(-4  2_q^^.^^^.g.g.)] 

=  (2Tr)"2^D"2|(h,  ,)h  ;dg,.../dg  s,g„...g  Z(g,h),  (2.12) 


-1  r  •   0 ""■'■   0^12  -  - "n 


where 


.-tn 


Z(g,h)=(2,r)-=-'|(h..)i-^exp(-i  I  q„^i,„^jgigj).  (2.13) 

1  »  J  —  i. 

Now  Z(g,h)  is  an  ordinary  normal  distribution  of  n  variables 
(giir''»Sj,)  with  the  c ©variance  matrix  (h.  •).  For  conve- 
nience, we  normalize  the  covariance  matrix  (h. .)  such  that 
the  (i,3)th  element  of  the  new  covariance  matrix  is 
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and  (v. .)  is  the  covariance  matrix  of  the  new  variables 
v/.=g./h''^.  Then,  with  the  above  change  of  variables, 
Equation  (2.12)  becomes  (Longuet-Hlggins) 

W(  +  ,,..,+)=(2Tr)-H-*|(h^.)|*(h^^h22...h^)*Jn,   (2.1^) 


where 


V  /dw3^...Xdw^w^W2...w^Z(w,v), 


(2.15) 


and  Z(w,v)  is  a  normal  probability  density  of  n  variables 
(v/^,...,w  )  with  the  normalized  covariance  matrix  (v.  •). 

By  Jac obi's  theorem  (Appendix  A),  the  determinant  and 
(i,j)th  element  of  the  matrix  (h. .)  are  computed  as 


j(h..)|=:D/D^, 


(2.16). 


_1 


^ir°i  ^ 


^11 '  ^In'  ^Ij' 


^nl ^rxn'  ^nj' 

-^n'''--"-^in'''^ij" 


(2.17) 


where 


°1= 


^11' ••••^In 


\l"-"^nn 


(2,18) 


Then,  by  substitution  of  Equation  (2.16)  into  Equation 
(2.1^), 
■       W(  +  ,,...,+)  =  (2Tr)-*"D^-^(h^^h22...h^)*J^.         (2,19) 
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We  nov/  compute  W  for  the  special  cases  when  n  is  1,  2 
or  3«  If  n=l,  then  Z(w,v);is  a  normal  distribution  of  a 
single  variate,  and  we  have 

J.=  /dw(2Tr)"2wexiD(-iw'^)=(2Tr)"2, 
^0 


d^=|RiJ=r(o). 


and 


hi,=D,-= 


^11'   ^^11 
-T?   '  -R   " 


=R(0) 


-1 


.R(0)  ,  R(0)' 
-R(0)',-R(0)" 


=-R(0)". 


Thus  Equation  (2.19)  for  the  case  n=l  becomes 


W(+)=(2Tf)-«D^-2h^^^J^ 

=  (2tt)"^[-R(0)"/R(0)]*. 


(2,20) 


When  n=2  or  3,  we  may  use  the  results  of  Nabeya(30) 
and  Kamat(19)  to  compute  J^^   and  J^: 


2A 


J2=(2-rr)""[(l-v^2^)2+v^2=os-^(-v^2)], 
J3=(2Tr)"-/^[|(v.^^)j*+(s^b^+S2b2+S3b^)], 


where 


-ir 


2^4 


,=cos-^L(v3,v,3-V33)(l-v3^2j-^(,_,^,^^2^-^-^^ 


(2.21) 
(2.22) 


^I=^3l''l2^^23- 


The  angles  of  arc  cosine  are  to  be  taken  in  the  range  (O.tt), 
and  Sg.  s^,  b^  and  b„  are  obtained  by  cyclic  permutation  of 


19 


the  V. ..  This  gives 

W(+.+)=(2^)-2(R^^R22-R^2^2^)-*(h^^h32)*[(l-v^2^)*4- 

•  v^^cos-^(-v^^)3,  (2.23) 

and 

W(+,+,+) 
=  (2Tr)~\"*(h^^h22h33)"^l|(v^j)|*+(s^b^+s2h2+S3b^)]. 

(2.2iJ') 
Since,  for  our  purpose,  we  need  to  compute  W(+,-)  and 
W(t,-,-),  we  consider  the  case  when  there  are  down-crossings 
in  W.  Suppose  the  kth  zero-crossing  is  down-crossing.  Then 
in  calculating  the  corresponding  probability  density  W  we 
need  only  to  take  the  range  of  integration  of  g^  from  -^   to 

0  instead  of  0  to  <».  Equivalently,  we  may  simply  reverse  the 

-1 
sign  of  the  (n+k)th  row  and  column  of  M  ,  and  hence  the  kth 

row  and  column  of  (h. .)  and  (v. .).  Thus  by  changing  sign  of 

X  J    -  J-  J 

v.g  in  Equation  (2.23)  we  obtain  (Longuet-Higgins) 

W(+,-) 
=  (2Tr)-2(R^^R22-R^2R3^)-*(h^^h22)[(l-vJ)*-v^2COS-lv^2^. 

(2.25) 

and  also  by  changing  corresponding  signs  of  v^.  in  s  and  b 
in  Equation  (2.23),  we  get 


W 


(,+,-, -)  =  (2TT)-^D^"^'h^^h22h33)^[  I  (v^^)|^+s^b^+ 

is^-Tr)h^Hs^-vYo^'].  (2.26) 
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Zero-crossing  Interval  Density  Function  for  a  Bandlimited 
White  Gaussian  Noise 


In  this  study  the  noise  is  assumed  to  be  ■bandlimited 
white  Gaussian  with  a  power  spectral  density 


S(w)=  j  -rr,   for  jwj<l, 
1  0,   for  |w|>l, 


(2.27) 


and  autocorrelation  function 


R(T)=r  sinr. 


(2.28) 


The  power  spectral  density  and  the  autocorrelation  function 
of  this  noise  have  the  shapes  as  shown  in  Figure  1 . 


S(w) 


-1         1 
(a) 


w 


(b) 


Figure  1. 

The  zero-crossing  interval  density  function  of  the 
noise  is  computed  according  to  Equation  (1.8) 

p(0;r)=W(+,-)/W(+)-   S      dt„[w(+,-,-)/w(+)] 
t^<t2<t3 

and  plotted  in  Figure  2.  Experimental  data  by  Blotekjaer(6) 
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•H 


22 

are  also  included  for  a  comparison.  The  programs  used  in  the 
computation  are  listed  under  LONGO,  LONGl,  and  MATIN.  For 
given  r I  W(+)  and  W(+,-)  are  computed  by  LONGO  according  to 
the  Equations  (2.20)  and  (2.25)  respectively  with  the 
autocorrelation  function  of  noise 

R(r)=-r~^sinT. 

The  program  LONGO  calls  the  subroutine  LONGl  to  compute  the 
term 

J-   dt^[W(+,-,-)/W(+)], 

where  T=to-t^. 

The  computation  repeats  with  different  value  of t  to  cover 
the  range  0<T<15.  MATIN  is  a  subroutine  called  by  LONGO  and 
LONGl  to  compute  the  determinants  of  matrices  which  are 
needed  in  evaluating  Equations  (2.17),  (2.18)  and  ((v. .)| 
in  Equation  (2.26). 
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Evaluation  of  W  for  Sine  Wave  Plus  Noise 

For  the  zero-crossing  interval  density  function  of  sine 
wave  plus  noise,  we  will  use  the  approximation 

p(0;r)=W(+,-)/W(+).  (3.1) 

The  assumption  used  here  is  that  the  signal-to-noise  ratio 
is  fairly  high  so  that  the  zero-crossing  intervals  are  mainly 
determined  by  the  signal  frequency,  and  the  probability  of 
having  two  down-crossings  in  the  interval t  is  negligible. 
Let  us  denote  the  signal  plus  noise  by  f(t).  Thus 

f(t)=s(t)+n(t),  (3.2) 

where  n(t)  is  stationary  Gaussian  random  noise  with  an 
autocorrelation  function  of 

R(r)=T  sinr, 

and  where 

s(t)=acos(qt+eQ),  with 

a=ratio  of  sine  wave  amplitude  to  rms  noise, 

q=radian  frequency  of  the  signal,  and 

eQ=initial  phase  angle  which  is  unknown. 
The  signal  s(t)  is  assumed  to  be  statistically  independent 
to  the  noise  n(t). 

To  compute  WC+)  we  consider,  as  before,  the  joint 
probability  that  f (t)  should  pass  zero  with  a  slope  g^  in 
the  time  interval  (t.,t.+dt.),  that  is,  the  probability  of 


zk 


f(1;)=s(t)+n(t)=0, 
f«<t)=s'(1;)+n'(t)=g^, 

or ,  equivalently , 

n(t)=-s(t), 


(t^<t<t^+dt^), 


(3.3) 


n«(t}=g3^-s'(t), 


(t^<t<t^+dt^). 


(3.^) 


Since  the  signal  and  noise  are  statistically  independent, 
the  proI)aMlity  that  the  conditions  in  Equation  (3 '3)  are 
satisfied  is  determined  "by  the  statistics  of  the  noise  alone; 
that  is,  it  is  equal  to  the  probability  of  the  noise  satis- 
fying tne  conditions  in  Equation  {3.'+).   Thus 


p(f^=0,f^'-g-j^)=p(n^=-s^,n^'=--g^-s^') 

=  (2Tr )  "^ I  M I  "^exp ( -Ix'^M'-'-X ) , 


(3.5) 


where 


x=r-s. 


-acosQ 


g. -s^ •  :    g. +qasin9 


M=fR^^. 


11 


i^%l*'"^ll' 


1,  0 
O.-Rq". 


9=the  instantaneous  phase  angle  of  the  sine  wave  given 

l)y-e-qt+6p, 
R=the  autocorrelation  function  of  the  noise. 

Since  6q  is  unknown,  we  may  with  a  sufficiently  large 
sampling  time  consider  9  as  a  random  variable  which  can  take 
all  values  from  -rr  to  rr  v,dth  equall  probability.  Hence  from 
Equation  (2.2) 


v'''^":^'....;"'2:5 

V/(-f-)=  (2tt 

-TT   0  , 

,') 

=  (2it) 

"^  i"deXdg.gj(2Tr)"^j 

-TT   0 

(- 

Mj  -^'expC-lCa^cos^e  + 
Ro")"^(Si+qasine)^]] 

-(2it) 

■^(-Rq")"*  /deexp(- 

2   2°° 
ia  cos  e.Udg.g. 

Q       1    i 

©xpC-K-Rq"  )~''"(gi-Kiasin 

e)2]. 

By 

a  change  of 

g=(-V>" 

variable 
^(gj_+ciasine) , 

W(+)=(2Tr) 

"^(-R  ")~*  Jdeexp(- 

-TT 

-i-a2cos2e)jrdg(-RQ")^[( 
c 

o 
qasin9]exp(-Jg  ) 

» 

1 

^g- 

where 

c=(-Ro"r^ 

^qasine. 

Letting  b=(-RQ' 

•)"*qa. 

<7 

W(+)  =  (2Tr) 

"3/2(_R^..)i;deexp( 

-ia^cos^e)[ 

(2Tr 

"^    2 

)  ^;   dggexpc-ig'^) 

bsine 

-bsine(2Tr)~2  S     dgexp 
bsine 

(-*g^ 

)| 

=  (2Tr)- 

^^  (-Ro")^Xdeexp(- 

|a^cos^e)[(2Tr)~*exp(- 

|-b^sin-6) 

-bsino[i-(2Trj~'^ 

bsine        -> 
S     dgexp(-#g2)]]. 
0 

1  X          p 

Since  erf(x)=2Tf~2  XdtexpC-t'^)  and 

n 

the  average  of  sin^: 

from.. 

-TT 

to 

TT  is  zero, 

\j 

W(T)  =  (2Tf)' 

■3/^(-R_")'^""  JdeexpC 

-TT 

-ia^cos^e)  [ 

.  '-  -.  .  , 

^ 

>.6 
(2TT)~^exp(-|b^sin^e)+ibsln6erf(2"2b3ine)].      (3.6) 

This  result  was  given  by  Rice (38). 

To  compute  W(+,-)  we  need  to  consider  the  instance  when 


or 


f(t)=s(t)+n(t)=0, 
f(t)=s»(t)+n'(t)=g, 


(t^<t<t^+dt^)  (i=l,2), 


n(t)=-s(t), 

(t.<t<t.+dt. )  (i-1,2). 
nMt)=g,-s'(t),  til   1 


(3.7) 


(3.8) 


The  probability  density  for  this  instance  is 


pCn^.ng.n^'.ng' )=(2^)" 


M 


■2exp(-|X^M~-'-X), 


(3.9) 


where 


X= 


-s. 

= 

-acose 

-^2 

-acosCe+qr) 

.^2"^2', 

g^+qasine 
g2+qasin(a+qr) 

M= 


^11'  ^12»  ^11'  ^12 

^21'  ^22'  ^21'  ^22 

p '   _p »   _p"   _P" 
~^11' ■12*  "11'  ^12 

P»      .l;P  '-     _P"      _P" 

""21 '  "22 '  ^-21'  *^22 


N 


y 


1  ,  R  ,  0  ,  R' 
R  ,  1  ,  R',  0 
0  .-R'.-R^.-R" 


-R',  0  ,-R"i-Ri; 


/ 


Oj 


Again  we  consider  0  as  a  random  variable  distributed 
evenly  from  -tr  to  tt.  And  recognizing  that  the  range  of  gp  is 
from  -00  to  0,  v/e  write  from  Equation  (2.2) 


Tf   ^"^ 


_1   "         " 

W(+,-)  =  (2Tf)   ..fdej'dg.j"  dg^g.  (-g,)p(n.  ,np,n  '  ,np' ) 
,.Tf   0    -oo   '^      ~     X   /C   X    ^ 


2? 


1     Tf       00  00 


where 


:(2Tf)"^|Mr^  Xd9;dgJ'dg„g.gpexp(-iuV^U),  (3.1.0) 


U= 


-acos9 
-acos(9+qr) 
g^+qasin9 
-g2+qasin(9+qr) 

The  exponent  in  Equation  (3.10)  may  be  written  as 

=-(2|M|  )"-'-^2a^M^^[cos^94-cos^(9+qr)]+2a^Mj^2j,'^osecos(9+qr) 
-2aM-.2[(gi +qasin9  )cos9+[g2-qasin(9+qr)]cos(9+qr)] 
+2aM. „[[g2-qaRin(0+qr)]cos9+(g.+qasin9)cos(9+qr)] 

: +M22C  (gi+cLasin9  )  +[g2-qasin.(9+qr)]~] 
-2M2o(g;j_+qasin9  )rg2-qasin(9+qr)]  y  , 


(3.11) 


where 


M^^-r (-Rq" )^-(-R" )^]-(-Ro" ) (R* )^  . 
M^2=R'[-(-Rq")R+(-R")], 

Ml3=-R'[(-Ro")-R(-R")]+(R')^, 


II  \2 


.^2 


M22=(-Ro")(i-R~)-(R')^, 

M23=-(-R")(i-R^)+R(R')^» 

!m|-[(i-r)[(-Rq")+(-R")]-(r')^][(i^r)1;(-Ro")-(-R")]-(R')^]. 


.\:.  28  ■  .•    ■      . 

We  make  a  change  of  variables 

(2|M|)'*g^=rcos(X+|Tf), 

(2lM|)"*g2=rsin(X+iti-), 
vmich  yields,  aftej?  some  rearrangement  of  terms, 

W(+,-)- 
=-iTf~^\M\^^'^   XdoXdr7dXr^ccs(2x)exp[-(D.r^+2D,r+D,J], 

-TT  ^  0  .  ..|lT  .  1       -     JJ 

(3.12) 
where 

D^=M22-M23COs(2x) 

=  (M22-M23)cos^\+(M22-HVl23)sin^X, 

D2=a|Mr*[-[(M^2-Ml3)'=os(iqr)+q(]\'[22-M23)sin(-|qT)] 

cos(9+i-qr)co3\ 
+[  (M^2"^^'^i3 )  ^^^^  s^"^^  ~^^^22'^^2  3  ^  °  °^  ^^^'''^  ]sin(  e+iqT)  sinx] , 

D^=a2|M|"^[[(M^^-H:^I^^)cos^(iq-r)+2q(M^2-^13)sin(iqT)cos(iqr) 
+q^  (M22-M23 )  sin^  (iqr)  Icos^  ( e+i-qr) 
+[(M^^-M^j^)sin^(iq-r)-2q(M^2+%3)sin(iqr)cos(iqr) 
+q^  (M22+M23  )cos^  (iqT)  ]sin^  ( 8+JqT)  ] . 

We  simplify  the  above  expressions  by  v/riting 

C^=l-R, 
C2=l+R, 
C3=-R^"+R". 


C^=(l-R) (-Rq"-R" )-(R' )^=C^C^-Cq^, 
C5=(l+R) (-Ro"+R" )-(R' )^=C2C3-Cq^, 
Cr,=CQC0s(iqT)+qC2sin(iqT) , 
Cg=CQsin(iqT)+q.C^cbs(iqr)  f 

and  replacing  (Q+iqr)  by  9  to  obtain 

D^=C2C^cos^\+C^C^sin^X, 

D2=~a(C^Cg)"^(C^Cr,cosecosX+C^Cgsinesinx) , 

D^=a^(C^C^)"^[C2"^C^(C^^4^gCOs^(iqT)cos^e 

+C^"-'-Cg(CQ^+C^sin^(Jqr))sin^9], 

|M|=C3C^. 

By  integrating  the  Equation  (3.12)  with  respect  to  r 
we  finally  write  (Cobb) 

W(+,-)=(2iT)'^|Mp/^  /cle  /dXD."^cos(2\)[(l+D,,^)exp(-D^) 

-TT^Dij,  (  3/2+D^2 )  ( 1 -er f Dj^ )  exp  ( -D^ )  ] , 

(3.13) 


D,=D,/D,S 


where 

V^2';"l 
D3=D3-D^2^ 


30 


Zero-crossing  Interval  Density  Function  for 
Sine  Wave  Plus  Noise 


The  zero-crossing  interval  density  function  for  sine 
wave  plus  noise  is  computed  according  to  the  Equation  (3.1) 

p(0;r)=W(+,-)/W(+) 
with  the  same  noise  characteristics  as  before.  Several  of 
these  density  functions  are  plotted  in  Figures  3-8.  STOSN, 
SINO,  and  COREK  are  the  programs  used  to  compute  p(0;t). 
The  computation  of  W(+,-)  is  carried  out  by  SINO  according 
to  the  Equation  (3- 13).  The  program  COREK  computes  Equation 
(3.6)  to  evaluate  W(+)  for  the  signal-plus-noise  case.  The 
data  computed  by  SINO  were  stored  in  the  computer  operating 
system  by  the  program  STOSN,  and  later  retrieved  by  COREK 
to  complete  the  computation. 

The  second  term  in  Equation  (1.8) 

;  dt  [W(+,-,-)/W(+)]  (3.1^) 

^1<V^3 

represents  the  probability  density  of  there  being  at  least 
one  down-crossing  between  an  up-crossing  and  a  down-crossing 
separated  by  an  interval  T=t^-t. .  When  the  amplitude  of  the 
signal  is  high,  the  zero-crossing  intervals  of  the  signal 
plus  noise  are  close  to  the  half  period  of  the  signal  wave, 
Tj^=Tr/q,  and  the  probability  that  the  zero-crossing  intervals 
are  much  longer  than  T,  is  very  low.  In  other  words,  if  t  is 
larger  thaji  T,  ,  the  probability  of  having  at  least  one  down- 
crossing  within  the  interval  '''  is  very  high. 

The  zero-crossing  interval  density  function  in  Equation 


.,■■.  31   ,-■  ■■^;. 
(3.1)  therefore  contains,  with  the  omission  of  the  second 
term  in  Equation  (1.8),  erroneous  distribution  near'r=3Tj^. 
When  the  signal  radian  frequency  q  is  0.6  or  higher,  the 
erroneous  distribution  fails  within  the  interval  (0,15)  of  r. 
These  erroneous  data  are  discarded,  and  the  main  distribution 
curves  are  extended  smoothly  to  vanish  at  r=15  by  the  program 
DPOL. 

The  zero-crossing  interval  density  function  (3.1)  is  a 
good  approximation  if  the  ratio  of  sine  wave  amplitude  to 
rms  noise  is  large.  In  the  cases  when  the  ratio  a=l  or  2,  the 
error  in  the  approximation  is  fairly  high.  Instead  of  compu- 
ting the  density  functions  for  these  cases  directly  from  the 
Equation  (3.1) »  we  interpolated  the  data  by  the  program 
SAPOL  from  the  density  functions  for  a=0,3,4  and  5  using  a 
spline  interpolating  function  (Appendix  D).  The  zero-crossing 
interval  density  function  for  a=0  is  just  the  density  function 
of  noise  alone. 

The  pr^ogram  INPOL  expands  the  data  points  of  the  density 
function  by  locally  interpolating  with  a  spline  function.  The 
area  under  the  computed  zero-crossing  interval  density  func- 
tion usually  has  a  slight  deviation  {±5fo)   from  1,  and  this 
is  corrected  by  normalization. 
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CHAPTER  ill 
BtPIfflyENTATION  AND  EVALUATION  OP  THE  DETECTOR 


Implementation  of  the  Detector 

In  impleiEenting  the  detector  it  is  assumed  that  the 
signal-to-noise  ratio  and  the  length  of  a  signal  are  known 
to  us.  We  also  assume  that  the  signal-to-noise  ratio  is  fairly 
large  so  that  the  additive  noise  on  a  signal  only  disturbs 
the  zero-crossing  interval  determined  by  the  signal  frequency 
and  does  not  add  more  zero-crossings  in  the  process.  Then  we 
know  the  number  of  zero-crossings  needed  to  observe  the  signal 
in  full  length. 

The  schematic  of  the  detector  is  shown  in  Figure  9. 
Zero-crossirigs  of  f(t),  both  up-crossings  and  down-crossings, 
are  first  detected  and  the  output  pulse  of  the  zero-crossing 
detector  is  used  to  strobe  the  system.  The  intervals  between 
zero-crossings  are  measured  and  compared  to  predetermined 
thresholds.  The  output  Z.  of  this  comparator  has  values  such 
as 

Z^=  {  1,        if  Tq<T^<T^, 

\  0,        otherwise,  (^.1) 

where 

TQ=the  lower  limit  of  the  threshold  of  the  zero-crossing 
interval,  comparator, 

38 


39 


•H 

o 

O-P 

N 

1   cs 

' 

PiO 

^^ 

•d 

. 

<D   O 
■P  .c! 

y 

r /^ 

fi  ra 

•H 

3   0) 

t<J 

o  y^ 

CJ 
1 

•P 

1 

•H 
CS3 

J^  ra 

0)   >s 

-p-p  oS 

V, 

. 

N. 

<H  ran 

V 

^  M-d 

•H 

W  0) 

Q 

5^  C 

cvj 

/ 

^ 

•H 

ta 

•H 

X 

Jh 

W) 

W)H  O 
1    S^csJ-P 

fi 

O  .H  >  (Jj 

•H 

Jh  ra  ^  f^ 

ra  J4 

0)  ra  <D  d 

CQ  O 

Nl  o  -P  ft 

O-P 

^4   fi  S 

U  O 

O  .H  O 

O  0) 

o 

1  -p 

O  0) 

/ 

V 

Q> 

N 

/ 

V 

■ 

•r 

> 

MH  ;.| 

^-N 

1   C  kJ  a> 

+> 

O'H  >  Jh 

^^ 

N. 

fj  ra  ^  5 

«H                                                            -^ 

c  to  o;  CQ 

^a  o-p  ctj 

f^  G  o 

O.H  s 

OS 


lai) 


T.=the  higher  limit  of  the  threshold  of  the  zero-crossing 
interval  comparator. 
"The  outputs  of  the  zero-crossing  interval  comparator  are 
stored  in  a  shift  register  and  the- content  of  the  shift  regis- 
ter is  counted  by  an  up-down  counter.  If  the  signal  s(t)  is 
composed  of  n  cycles  of  sinusoidal  waves,  then  there  are  2n 
zero-crossings  generated  by  the  signal.  To  observe  the  signal 
in  full  length,  the  counter  sums  up  the  content  of  the  shift 
register  in  such  a  way  that 

'0,         if  i=G, 

if  z.=z;: 


C^_^-l,    if  Z^=0  and  Z^_2^=l,. 


where  Z.=0  for  i^O.  Therefore  the  value  of  C.  at  any  moment  is 

C.=  2   Z,  (if. 3) 

^   k=i-2n  ^ 


and  this  is  compared  to  a  threshold  to  determine  the  presence 
of  a  signal.  Thus  if  we  set  the  threshold  to  m,  then  the 
output  of  the  comparator  is 


D^=  I  1,   if  m^C^^2n, 

\  0,   otherwise.  i^A) 


Let  us  denote  that 

H^=the  hypothesis  that  a  signal  is  present, 
HQ=the  hypothesis  that  a  signal  is  absent, 
and  for  given  a  and  q  of  the  signal, 


^1 

P  =?(Z.=l|H^ )=the  probability  that  when  a  signal  is 

present  T.  lies  within  the  interval  (Tq,T. ), 
Q  =P(Z.=0|H. )=the  probability  that  when  a  signal  is 

S      J-      JL 

present  T.  lies  outside  the  interval  (T„,T^), 
P  =P(Z.=l|HQ)=the  probability  that  when  a  signal  is 

absent  T.  lies  within  the  interval  (Tq,T^), 
Q  =P(Z.=0|H„)=the  probability  that  when  a  signal  is 

absent  T.  lies  outside  the  interval  (Tq,T^). 


Then 


and  P  corresponds  to  the  area  under  the  curve  bounded  by 
T=T^  and  T.  in  Figure  10(b),  and  P  corresponds  to  the  area 
under  the  curve  bounded  by  the  same  limits  in  Figure  10(a). 
To  simplify  the  analysis,  the  detection  probability  of 
a  signal  is  approximated  by  the  detection  probability  deter- 
mined by  comparing  threshold  at  the  end  of  the  signal..  This  is 
the  lower  bound  of  the  detection  probability,  and  is  given  by 

P^=P(D^=:l|H^)=P(C^^m|H^) 

k=m 

The  corresponding  false  alarm  rate  is 

Pf=P(D.=l|HQ)=P(C.^m|HQ) 

2n 

k=m  "^  "  " 


^2 


<a>- 


Cb) 


Pigiire  10, 


^3 

Now  the  objective  is  to  maximize  the  detection  proba- 
bility P,  for  a  given  false  alarm  probability  P^.  This 
strategy  is  known  as  the  Neyman-Pear s on  criterion  of  decision 
making,  which  does  not  involve  either  a  priori  probability 
or  cost  estimates  (Whalen) . 

Let  r=P  /P^.  Then 


s'^  n 


2n 
k=m 


The  binomial  expansion 


n 

L 

k=ra 


f(n,m;x)=  I  (5J)x^(l-x)^-^  (^.8) 


is  related  to  the  incomplete  beta,  function  by 

f(n.m;x)=r  /dtt"^-^(l-t)^-"^]/C  )dtt^-l(l-t)^-"^J    (^'9) 
0  0 

(Abramowitz) .  Since  the  right  side  of  the  above  equation 
increases  raonotonically  with  respect  to  x  for  given  values 
of  n  and  m,  the  binomial  expansion  on  the  left  side  also 
monctonically  increases  with  respect  to  x.  Then  we  see  that, 
for  a  signal  with  parameters  a,  q  and  n,  P^  in  Equation  (^.7) 
is  maximized  by  maximizing  the  ratio  r  for  given  m  and  P^, 
and  thus  P^..  The  maximization  of  r  is  done  by  the  program 

rcoM3. 

PC0M3  integrates  the  zero-crossing  interval  density 
function  of  noise  from  r-0  until  the  integrated  area  equals 
P  ,  and  the  zero-crossing  interval  density  function  of  signal 


plus  noise  is  also  integrated  for  the  same  interval  to  obtain 

the  corresponding  P  .  This  integration  procedure  is  repeated 

with  an  incremented  lower  limit  of  integration  until  the 

range  of  r  is  exhausted.  Among  the  P„'s  obtained  by  such 

s 

integrations,  the  largest  value  is  picked  out  as  the  maxi- 
mized P  for  a  given  P  .  The  integration  interval  correspond- 
s  n 

ing  to  the  maximized  P  is  then  the  optimum  threshold  for 

s 

zero-crossing  interval  as  in  Equation  (^.1). 

Figures  11  through  1^  show  the  relations  between  P  and 
the  maximized  P  for  different  values  of  parameters.  In  the 
Figiires  13  and  1^,  P  is  plotted  against  q  with  fixed  para- 
meters P  and  a.  From  these  graphs  we  observe  that  q=0.4  in 

general  gives  a  large  value  to  the  ratio  r=P  /P  . 

s  n 

The  final  step  of  maximizing  the  detection  probability 
P^  for  given  P^  is  the  determination  of  the  threshold  m.  The 
program  RCOM  computes  P^  and  P^  for  different  values  of  P 
and  m  according  to  the  Equations  (^.7)  and  (4.6)  respectively 
as  shown  in  Figures  15-19.  From  these  graphs  the  optimum 
value  for  m  could  be  determined  to  give  the  maximum  P,  for 
given  P^  of  a  signal  with  given  parameters.  For  instance, 
the  optimum  value  of  m  for  a  signal  with  parameters  a=2, 
q=0.4  and  n=4  to  be  detected  with  maximum  probability  for 
a  given  false  alarm  probability  P^=10~^  is  6,  while  7  is  the 
optimum  for  P„=10  . 
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Comparison  to  Correlation  Detector 

The  performance  of  the  detector  is  compared  to  the 
conventional  correlation  detector.  Suppose  that  the  para- 
meters of  the  signal  s(t)  are  completely  known  to  us,  and 
let  H.  and  H„  denote  the  hypotheses  that  the  signal  is 
present  and  absent  respectively.  The  correlation  detector 
chooses  H.  if 


/f(t)s(t)dt^V, 
0 


(5.1) 


where  the  interval  (0,T)  is  the  duration  in  which  the  signals 
are  received  and  V  is  a  threshold.  The  implementation  of  the 
correlation  detector  is  shown  in  Figure  20. 


'0  or 


Figure  20. 


A  correlation  detector  is  equivalent  to  a  matched  filter 
followed  by  a  threshold  comparator  (Whalen) .  This  equivalent 
detector  which  implements  the  operation  of  the  correlation 


55 
detector  using  a  matched  filter  with  impulse  response  h(t) 
is  shown  in  Figure  21. 


f(t) 
O^t^T 


h(t)=s(T-t) 
0<t^T 


^ 


Sampled  at 
t=T 


Threshold 


Hq  or 


Figure  21 . 


The  decision  rule  in  Equation  (5.1)  stands  on  the 
assumption  that  the  noise  is  white  Gaussian.  ..The  -threshold 
V  is  determined  to  meet  an  employed  detection  criterion  for 
the  particular  application.  To.mak^-  the  environment  similar 
to  our  previous  analysis  in  analyzing  the  correlation  detec- 
tor, we  assume  that  the  noise  is  handlimited  white  Gaussian 
with  power  spectral  density 


S(w)=  J  Tf,   for  |w|^l, 
0,   otherwise 


(5.2) 


and  autocorrelation  function 


R(T)=--r  sinr. 


(5.3) 


The  Neyman-Pearson  criterion,  which  involves  neither  a  priori 
probabilities  nor  cost  estimates,  is  also  employed  here.;  The 
strategy  of  this  criterion  is  to  m.aximize  the  detection 
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probability  for  a  given  false  alarm  probability,  and  this 
could  be  accomplished  by  using  a  likelihood  ratio  test 
(Whalen).  We  approach  the  analysis  with  the  case  where  the 
signal  is  sampled  at  m  discrete  times. 

Let  us  denote  the  ra  sampled  values  by 

fj^=S3^j^+nj^,  (k=l,.  ..,m)  (i=0,l),  (5.^) 

where 

^ok=°'  ^ik=^(\^'  V^^V'  V^^\^' 

and  t,  (k=l,...,m)  are  sampling  times. 

If  a  noise  is  sampled  at  in  discrete  points,  the 
covariance  matrix  M=(m.  .)  of  the  noise  has  the  elem.ents 

m^^=E[n^np=R(t^-t^.)  (i,j=l m).  {5.5) 

Since  the  noise  is  Gaussian,  the  joint  probability  density 
of  the  m  sampled  values  of  noise  is 


p(n)=p(n^,n2,...,n^) 

=  (2tt)~'^''"|M|  -*exp(-|xV^X),  {5.6) 

where  X  is  the  column  vector  of  sampled  values.  By  writing 

Q=M"^=i((1^^0Mi,j-l m), 

the  Equation  {5.6)   becomes 

p(n)=(2TT)"^"'|Mr"2'exp(-|  S   q.  .n.n.).  (5.?) 

i,D=l   ''  ^ 

Since  n,=f,-s,,  (h--0,l),  the  above  equation  could  be 


51 
vnritten  as 

p(n)=p(f-s^) 

=  (2-rr)-*"^|M|-*exp[-i_    2^  ^ij^^i-^hi^  ^^j-^hj^^  V 

Then,  when  the  signal  is  absent  in  f (t) ,  h=0  and  the  joint 
probability  density  of  m  sampled  values  o^  (^^"^ok^ 
(k=?l, . .  .  ,m)  is 

PQ(f)=p(f-SQ)=p(n) 

A     i       m 

(5.8) 

If  the  signal  is  present  in  f (t),  then  h=l  and  the 
corresponding  joint  probability  density  of  the  m  samples 
(fjj-s^^)  (k=l,  .  .  .  ,m)  is  , 

p^(f)=p(f-s^)=p(n) 


{5.^) 

The  likelihood  function  for  the  hypothesis  Hq  is 
the  joint  probability  density  function  p^Cf )=pQ(f. , . . . ,f^) . 
Likewise  p. (f^,...,f  )  represents  the  likelihood  function 
for  the  hypothesis  K^.  In  terms  of  these  likelihood  functions, 
a  likelihood  ratio  r  is  defined  by 

^=Pl(^l ^m^/Po^^l =^m) 

m 
=exp[-i.  S^^q,.(f.-s^,)(f.-s^.)3/ 


m 
E 
i,j=l  "^ 


exp[-i^  E_^q..(f.-So.)(f.-So^)1 
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m 
=exp[-i-^  E  q^-.(2f.So.-2f.s^.-So.So.4-s^.s^.)].    (5.10) 

and  the  likelihood  ratio  test  is  to  choose 

H^,  if  r^r^, 

Hq,   if  r<rQ,  (5.11) 

v/here  r^  is  the  decision  threshold. 

If  the  noise  is  not  white,  the  evaluation  of  the  likeli- 
hood ratio  r  in  Equation  (5.10)  involves  a  matrix  inversion 
from  the  covariance  matrix  M.  When  the  signal  is  sampled  at 
a  large  number  of  points,  this  matrix  inversion  is  not 
practical  in  considerations  of  speed  and  simplicity  of  the 
system.  This  is  also  the  case  with  our  bandlimited  white 
noise,  and  we  simplify  the  analysis  by  sampling  f(t)  at  the 
time  intervals  for  which  the  sampled  noises  are  uncorrelated. 

Since  the  autocorrelation  function  of  the  noise  has 
zeros  at  T=k.'rr  (k=±l,±2,  .  . . ) ,  if  ws  sample  f(t)  at  the  inter- 
vals 6t=Tr,  the  samples  are  uncorrelated;  that  is 

E[n^np=R(t^-tj^)=  [  R(0),   if  i=j 

I  0,      otherwise.  (5.12) 

Then,  since,  from  Equation  (5.3).  R(0)=1,  the  elements  of 
the  matrix  M  in  Equation  (5 '5)   become 

in..=E[n.n.>  ll.   if  i=3. 

I  0,    otherv/ise,  (5.13) 

and  thus, 

M=I, 
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and 

q..=  f  1.    if  i=j. 

I  0,   otherwise.  (5'1^) 

Then  the  likelihood  functions  in  Equation  (5.8)  and 

(5.9)  simply  become 

i„       m        5 
PQ(f)  =  (2iTr"^"™exp[-i  L  (-V^Ok^^  ^^-^-^ 

and 

1       m        ^ 
p^(f)=(2Tr)-^%xp[-4.L  (fk-Sii,)'^].  (5.16) 

The  likelihood  ratio  r  in  this  case  is 

r(f)=p^(f)/pQ(f) 

m        Q        m    .0 
k=l  ^  "^         k=l  ^  ^^ 
or,  upon  combining  terms, 

r(f)=exp[-i^L^(2Voj^-2f^s^^-So/+Si/)J.       (5.18) 

The  decision  rule  is  to  choose  H^^  if  ^^^Qt    or  equivalently, 
by  taking  natural  logarithm,  choose  H^  if 
m  99 

j^^,^Vik-:^k^ok*ok  -*^ik  >^i"(^o)-  (5-19) 

Let 

m  00 

%f/fk^lk-V0k^^0k  -*^lk  )■  (5-20) 

To  evaluate  the  performance  of  this  detector  we  need  to  know 
the  distribution  of  G  under  each  hypothesis.  When  Hq  is  true, 
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f  =nj^  and  the  expectation  of  G  is 


Eo[G]=E[  E    in^s^^-n^s^^Hs^^'^-hs^^^n.  (5.21) 

K — 1 


Since  Sqj^=0, 


m 


111  r\ 

Eo[G]=-i  Z  s^/.  (5.22) 

xC — 1 

The  variance  of  G  under  Hq  is 

Vor.G>E[(G-Eo[G])2].  (5.23) 

But 

m 

so 

mm 


VorG>E[  Z       E  n^n.s.^s^.] 
"  k=l   ,1=1       .'^  ""  ^ 


m       m 


Since 


E[nj,np='    !  R(0)=1,        if  k=j, 

I  0,  otherwise, 

it  becomes 
Define 

m         r) 

S=  Z   s,/.  (5.25) 

k=l  ^^■ 

Then, 

EqCg]— is, 


61 


Vq[g]=S. 


Therefore  the  density  function  of  G  is 

PQ(G)  =  (2Tr)~*S-*exp[-i(G+iS)Vs].  (5.26) 

In  a  similar  fashion  we  obtain  the  density  function  of 
G  under  the  hypothesis  H.  .  If  the  hypothesis  H.  is  true, 
then  ^k=s^, +n,  and 

m  m 


E^[G]=E[^S^(s^j^4-nj^)s^^+i^E^(-s,j^'^)] 


m    p   m    p 
k=l  ^^   k=l  ^^ 

=i  2  s.,'^.  (5.27) 

k=l  ^^ 

m  m 


i-E^[G]=^L^(s^^4«^)s,j^+i^E^(-s,^^)-4^E^E=^^' 


k=l  ^  ^^ 


V^[G]=E[(G-E^[G])^] 


m   m 

i:   L 

k=l  j=l 

m. .,  m 
J 


=^t;,.2;  Z\^k"j^ik^iP 


\?1  .^l^K^p^lk^lj 


m    5 
=  2  s.,-^.  (5.28) 

k=l  ^^ 

Thus 

p,(G)  =  (2Tr)"*S~WL-i(G-|-S)Vs].  (5-29) 
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Since  we  choose  H.  if  G^ln(rQ),  by  writing  c=ln(rQ), 
the  false  alarm  probability  is 

00 

P^=P(D^JHq)=  /dGpQ(G) 


1   1 


2 


=  (2tt)~2S"2  /dGexp[-i(G+iS)VS.]. 


c 


By  a  change  of  variable  z=S~^(G+iS), 


i  i   ~     i     1  2. 


P^=(2tt)"^S'=   J"  dzS^exp(-|z'=') 


(c+|S)/S^ 

=  (2Tr)"'^   /  ^dzexp(-iz^).  (5.30) 

(c+is)/S* 

The  probability  of  detection  is 

I 

00  ■' 

P^=P(D^|H^)=  J"dGp^(G) 

c  . 

=:(2'rr)''^S"^  XdGexp[-i(G-4S)Vs] 
c 

=  (2tt)"^S~^   J*  ^dzS^exp(-iz'=^) 
(c-|-S)/S* 

=  (27t)'2   X  ^dzexpC-iz'^)  (5.31) 

(c-iS)/S* 

The  performance  of  the  correlation  detector  is  computed 
by  the  program  PLIN  and  shov/n  in  Figure  22.  We  note  that  the 
performance  of  the  correlation  detector  is  superior  to  the 
zero-crossing  measuring  detector  in  the  assumed  circumstances 
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of  the  signal  and  background  noise.  For  instance  when  a=2, 
q=0.4  and  n=^,  the  false  alarm  probability  of  the  correla- 
tion detector  at  P,=0.9  is  ^xlO"'  while  the  false  alarm 
d 

probability  of  the  zero-crossing  measuring  detector  at. the 
same  detection  rate  is  3x10  '^ . 

This  superior  performance  of  the  correlation  detector 
to  the  zero-crossing  measuring  detector  is  based  on  the 
following  assumptions- about  the  signal  and  noise;  the 
frequency  and  amplitude  of  the  signal  are  completely  known 
and  fixed;  the  noise  is  stationary  in  Gaussian  distribution. 
If  the  practical  situation  deviates  from  the  assumptions, 
then  the  performances  of  the  two  detectors  will  also  deviate 
from  the  previous  results.  The  merits  and  operating  charac- 
teristics of  the  two  detectors  are  further  discussed  in  the 
next  chapter  under  such  practical  considerations  in  detecting 
signals  from  sleep  electroencephalograms (EEGs ) . 
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CHAPTER  IV 
-FURTHER  DISCUSSIONS  AND  CONCLUSION 


Operating  Characteristics  in  Practical  Situation 

The  zero-crossing  measuring  detector  has  been  applied 
to  the  detection  of  phasic  events  in  sleep  electroencephalo- 
grams(EEGs).  Examples  of  EEC  activities  from  a  cat  are  shown 
in  Figure  23.  The  arrows  indicate  phasic  events  in  a  sleep 
EEG.  Among  various  phasic  events,  alpha,  sigma,  and  beta 
waves  could  be  modeled  as  signals  of  sine  wave  form  imbedded 
in  the  noise  of  background  EEG. 

The  physiological  origins  of  the  phasic  events  as  well 
as  the  irregular  activities  of  EEG  have  not  been  completely 
understood  as  yet  and  neither  the  statistical  characteristics 
between  the  phasic  events  nor  the  irregular  activities  have 
been  studied.  Even  if  we  can  assume  for  practical  purposes 
that  the  phasic  events  are  independent  processes  superimposed 
on  the  background  EEG,  the  nonstationarity  and  the  lack  of 
Gaussian  property  in  the  sleep  EEG  activities  make  it  diffi- 
cult to  apply  the  previous  analysis  to  the  quantitative 
determination  of  the  advantages  of  the  two  methods  in  detect- 
ing phasic  events  from  sleep  EEG.  Hence  in  the  following  we 
attempt  only  to  make  qualitative  observations  on  the  merits 
of  the  zero-crossing  measuring  detector. 

(>5 


66 


CM 


P4 


','..:    67 

When  we  attempt  to  detect  phasic  events  from  a  sleep 
EEG  we  encounter  the  the  following  situations:  1).  The 
frequency  of  the  phasic  events  varies  among  individuals.  The 
frequency  of  alpha  activity,  for  instance,  ranges  between  8 
and  12  Hz,  sigma  activity  between  12  and  15,  and  the  a  priori 
probability  density  of  the  frequency  distribution  is  unknown. 
2).  The  amplitude  of  phasic  events  varies  due  to  the  elec- 
trode resistance  and  individual  differences  in  age  and 
physical  parameters.  3).  Background  EEG  activity  level  varies 
and  artifacts  are  introduced  by  movements.  The  background 
EEG  activity  varies  also  due  to  the  differences  between 
individuals  and,  within  an  individual,  the  variation  of 
activity  inherent  to  a  sleep  process.  Examples  of  movement 
artifacts  recorded  from  a  cat  are  shovm  in  Figure  26. 

To  qualitatively  observe  the  operating  characteristics 
of  the  two  detectors  in  the  described  situations,  we  do.  the 
following  analyses: 

1).  To  see  the  effect  of  the  varied  frequency  of  the 
signal  on  the  performance,  we  compute  the  detection  proba- 
bilities of  the  detectors  to  signals  of  different  frequency. 
Suppose  the  correlation  detector  is  tuned  to  signal  s-,  and 
we  receive  a  signal  Sp  With  a  same  amplitude  but  different 
frequency  from  s^.  Then,  since  ^y'^z'k^^'k.*   ^  ^^  Equation 
(5.20)  becomes 

m  m    p 

^\~/'2k-^''k)'lk-\^i"lk  •  ^^'^^ 

The  expectation  of  G  is 


^2^'^>J,^2tHk-iJ.Hk-  <6-2) 
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m         m 
k=l  "^^  ^^  k=l 
The  variance  of  G  is 

m   m 
k=l  D=l 

=  2  s.  J^=S.  (6.3) 

k=l  ^^ 

Thus  the  density  function  for  G  is 

P2(G)=(2Tr)"*S-*exp[-i(G-E2)Vs],  (6  A) 

where  E2=E2[g]. 

Since  the  detector  is  tuned  to  the  signal  s^,  the 
threshold  remains  the  same  as  for  s^ ,  and,  hence,  the  detect- 
ion probability  for  Sg  is 


P^=P(DjH2)=  /dGp2(G) 
c 


=  (2Tf)"2s"2\fdGexp[-i(G-E„)VS] 
c  . 

=  (2Tr)~*  S   ^dzexp(-iz^),  (6.5) 

(c-EgVS* 

where  H^  is  the  hypothesis  that  signal  S2  is  present,  and  c 
is  the  detector  threshold  set  for  signal  s. . 

The  detection  probability  of  the  correlation  detector 
for  signals  with  a  different  frequency  is  computed  by  the 
program  PSEG  with  the  following  parameters: 

a=2  for  both  s^  and  S2, 

n=6. 
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c=threshold  set  to  detect  s.  with  a  probability  of  0,95, 
that  is,  set  for  P(D^|  H^ )=0.95, 
where  a  is  the  ratio  of  signal  amplitude  to  rms  noise,  and 
n  is,  the  length  of  signal  in  niunber  of  cycles.  The  result  is 
shown  in  Figure  2^ (a).  The  detection  probability  for  a  signal 
with  a  different  frequency,  P(D^\E^),    sharply  drops  off  as 
the  frequency  of  the  signal  slightly  departs  from  the  tuned 
frequency  of  the  detector. 

The  detection  probability  PCD^lHg)  for  the  zero-crossing 
measuring  detector  is  computed  in  the  following  way.  Since 
the  detector  is  tuned  to  s^,    the  thresholds  m  in  Equation 
i^A)   and  Tq  and  T^  in  Equation  (^.1)  are  already  determined. 
From  the  Figures  16(a)  and  (b),  we  note  that  the  optimum 
threshold  for  m  is  6  for  a  signal  with  a=2  and  n=k,   and  m=12 
when  a=2  and  n=8.  So  with  the  present  signal  parameters  of 
a=2  and  n=6,  we  assume  that  the  threshold  for  m  has  been 
determined  at  9.  As  before  we  want  P(D^[  H^)=0.95,  which 
consequently  determines  the  optimim  thresholds  (Tq,T^)  for  the 
zero-crossing  interval. 

The  probability  of  detecting  s^   is  then 

Pd=P(Di|H2) 

2n  o 

\^„<  S'^z'd-^z)'"-'.  (6.6) 

where  Pg  is  the  area  under  the  zero-crossing  interval  density 
curve  of  Sg  plus  noise  bounded  by  the  interval  (Tq,T.,).  The 
computation  is  carried  out  by  the  program  PCOM^  and  RCOM, 
and  the  result  is  plotted  in  Figure  24(b). 


70 


1     : 
1 

i 

I 

« 

«' 

I 
7 

1 

i 

i 

N 

1 

1 

o 

o 

'o 

'o 

'o 

O 

'o 

'o 

'o 

'o 

.»_ 

11 

•  o 

A) 

M 

oP 

:  O 

N     <D 

"tJ 

:  0) 

II    R 

D    G 

U 

'.2 

•    .  71 

We  note  from  the  two  curves  that  the  operating  charac- 
teristic of  the  zero-crossing  measuring  detector  is  more 
suitable  to  cover  the  broad  range  of  signal  frequencies.  To 
cover  the  same  range  with  the  correlation  detector  without 
affecting  the  detection  and  false  alarm  performance,  we  need 
several  channels  of  detectors  in  parallel,  which  is  not  very 
practical. 

2).  The  zero-crossing  interval  distribution  is  not 
dependent  upon  the  amplitude  of  the  signal  and  noise  so  long 
as  the  signal  to  noise  power  ratio  is  maintained  constant, 
which  is  the  usual  case  with  sleep  EEGs.  And,  therefore,  the 
operating  characteristic  of  the  zero-crossing  measuring 
detector  is  not  affected  by  the  varied  strength  of  signal 
and  noise. 

However,  a  reduced  amplitude  of  a  signal  greatly  decrea- 
ses the  detection  probability  of  the  correlation  detector. 
The  detection  probability  of  the  correlation  detector  for  a 
reduced  signal  can  be  computed  in  the  same  way  as  done  for 
varied  frequency.  In  this  case  the  first  term  in  Equation 
(6,2)  becomes  a  function  of  amplitudes  of  two  signals  instead 
of  frequencies,  and  the  detection  probability  for  the  reduced 
signal  is  expressed  again  by  the  Equation  (6,5)   which  could 
be  computed  by  the  program  PSEG  with  a  few  modifications. 

The  changes  in  performance  of  the  two  detectors  are 
shown  in  Figure  25(a).  The  horizontal  coordinate  is  the  ampli- 
tude ratio  of  the  reduced  signal  s^   to  the  tuned  signal  s.  . 
The  frequency  and  the  length  of  Sp  are  assumed  to  be  the  same 
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as  s^  and  the  signal-to-noise  ratio  is  assumed  to  be  constant, 
The  detection  probability  of  the  correlation  detector  for 
signals  of  reduced  amplitude  drops  off  rapidly  as  the  ratio 
decreases  while   that    of  the  zero-crossing  measuring 
detector  remains  constant. 

3).  The  increased  noise  level,  assuming  that  the  noise 
characteristics  remain  the  same,  does  not  increase  the  false 
alarm  rate  of  the  zero-crossing  measuring  detector.  Also, 
since  the  signal  to  noise  power  ratio  usually  remains 
constant,  the  detection  probability  of  the  signal  is  not 
affected. 

However,  if  the  noise  level  is  increased  in  the  correl- 
ation detector,  the  false  alarm  rate  also  increases  to  a 

considerable  degree.  Suppose  the  power  level  of  the  noise 

2 

has  been  changed  by  a  factor  of  h  .  Then 

E[nj^n^]=  f  h^,   if  k=j 

\^0,    otherwise,  (6.7) 

and,  since  fj,=ni^. 

m    « 
E[G]=-i  E  s,^^=-iS,  (6.8) 


k=l  !>= 


P     mm 


V[G]=E[(G-E[G])^]=E[^_2^  ^E^nj^n^s^j^s^^] 


k=l  j=l 

2  °^    2   2 
=h^  E  s./=h^S.  (6.9) 

k=l  ^^ 
Thus 

PQ(G)=(2iT)"V^S"*exp[-i(G+|S)V(h^S)],  (6.10) 


7k 
and  the  false  alarm  probability  is 

00 

P^=P(D^JHq)=  /dGpQ(G) 
c 

=  (2TT)"V^S"*-XdGexp[-i(G+iS)V(h^S)] 

=  (2Tr)~2   X     dzexpC-iz'^),  (6.11) 

(c+iS)/(hS*) 

where  c  is  the  threshold  set  for  the  desired  detection 
probability  of  the  signal.  The  effect  of  the  increased  noise 
level  on  the  false  alarm  probability  of  the  correlation 
detector  is  computed  by  the  program  PNOI  according  to  the 
Equation  (6.11),  and  plotted  in  Figure  25(b).  The  false  alarm 
probability  of  the  correlation  detector  increases  rapidly  to 
0.1  and  gradually  approaches  to  0.5  as  the  noise  level 
increases.  The  false  alarm  probability  of  the  zero-crossing 
measuring  detector  remains  constant  regardless  of  the  increas- 
ed noise  level. 

Since  the  false  alarm  probability  of  the  correlation 
detector  is  susceptible  to  the  noise  level,  the  detection 
threshold  of  the  correlation  detector  could  not  be  lowered 
arbitrarily  to  increase  the  detection  probability  of  weak 
signals.  Consequently,  the  optimum  threshold  of  the  correla- 
tion detector  for  detecting  phasic  events  in  EEG  should  be 
detemiined  according  to  the  individual  strength  of  the  EEG 
activity. 

On  the  other  hand,  the  zero-crossing  measuring  detector 
does  not  need  any  threshold  adjustments  even  if  the  strength 
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of  the  EEG  activity  varies  over  a  wide  range  due.  to  indi- 
vidual differences  or  electrode  resistance. 

Typical  examples  of  the  detection  performance  of  the 
zero-crossing  measuring  detector  are  shown  in  Figure  27.  A 
square  pulse  is  the  output  of  the  detector  indicating  the 
presence  of  a  signal.  The  immunity  of  the  detector  to  high 
power  artifacts  is  demonstrated  in  Figure  26. 

Conclusion 

The  performance  of  the  zero-crossing  measuring  detector 
in  an  idealized  environment  of  signal  and  noise  was  inferior 
to  that  of  the  correlation  detector.  But  when  the  conditions 
of  the  signal  and  noise  deviated  from  the  assumptions,  the 
zero-crossing  measuring  detector  maintained  the  quality  of 
performance  while  the  performance  of  the  correlation  detector 
was  susceptible  to  the  change  of  environment.  The  sleep  EEG 
is  a  very  complex  process;  both  the  signal  and  noise  have  a 
wide  range  of  variation.  When  the  zero-crossing  measuring 
technique  is  employed  to  detect  phasic  events  from  sleep 
EEGs,  it  should  provide  a  detector  which  is  excellent  in 
performance  and  applicable  to  a  broad  range  of  subjects  with 
a  simplicity  of  implementation. 
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APPENDICES 


APPENDIX  A 
JACOBI'S  THEOREM 


In  evaluating  h. .  in  Equation  (2.17)  we  used  Jacobi's 
theorem  which  is  stated  as 

(A-^)^^^=|A|-^adj^^V.  (7.1) 

The  derivation  of  the  theorem  is  quite  lengthy,  and  so  we 
will  add  here  only  a  brief  explanation  on  the  notions  in 
Equation  (7-1)  and  the  application  of  the  theorem  in  evalu- 
ating h.  ..  A  reader  interested  in  the  development  of  the 

J 

theorem  could  refer  to  Aitken(2), 

Let  A  he  a  matrix  of  order  nxn.  The  determinant  obtained 
by  suppressing  m  rows  and  m  columns  of  jA| is  called  a  minor 

of  |Al  of  order  n-m. 

fk) 

The  kth  compound  A^  '   of  A  is  obtained  by  forming  its 

elements  with  minors  of  |a|  of  order  k  in  the  following  way; 
let  all  minors  which  come  from  the  same  group  of  k  rows  (or 
columns)  of  A  be  placed  in  the  same  row  (or  column)  of  this 
derived  matrix,  and  let  the  priority  of  elements  in  rows  or 
coliimns  of  this  matrix  be  decided  on  the  principle  by  which 
words  are  ordered  in  a  dictionary  or  lexicon.  For  example, 
minors  from  rows  1,2, A-  of  A  will  appear  in  earlier  rows  than 
those  from  1,2,5  or  l,3i^  or  2,3,^;  and  similarly  for  columns. 
Consider  a  minor  jB|  taken  from  rows  i.,i2f....i  and 
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columns  j^,J2»'''»J  of  k.    Then  the  complementary  minor 
formed  by  the  remaining  n-m  rows  and  n-m  columns  of  A,  with 
the  sign  factor  (-D^l"*"- •  •■^^m+Jl+' •  ■+Jm  multiplied  in  front, 
is  called  a  cofactor  of  minor  |Bl  . 

Then  the  kth  ad jugate  compound  of  A  which  is  denoted  hy 

(k) 
adj^  'A  is  obtained  in  the  following  way:  take  the  kth  com- 

(Ir) 

pound  A^   ,  replace  every  element  in  it  by  its  cofactor  in 
|A| ,  and  transpose  the  resulting  matrix. 

Now  let  the  matrix  M  in  Equation  (2.8)  be  equal  to  A 
in  Equation  (7.1).  Then 


-1 


M(^)=|M-V^adj(kV^ 


(7.2) 


and  the  determinant 


^1= 


^11 ^In 


^nl"-"^nn 


(7.3) 


is  the  element  in  the  1st  row  and  1st  column  of  the  nth 
compound  matrix  of  M.  This  element,  if  divided  by  |m|,  should 

be'  by  the  Jacobi   theorem  equal  to  the  element  in  the  1st 

■J 
row  and  1st  column  of  the  nth  ad jugate  compound  of  M~'  which 

is 


V 


"^n+l,n+l"-""^n+l,2n 


%,n+l'--"^2n,2n 


(7.^) 


where  the  matrix  (m.  .)=M~  .  Thus 

J.  J 


81 


D^=\iii\D^=m^. 


Since 


|(h,.)|=D2-S 


(7.5) 


(7.6) 


by  substituting  Equation  {7-5)   into  Equation  (7.6)  we  obtain 
the  relation  in  Equation  (2.16) 

Now  let 


P=(p.  .)  = 


n+l,n+l" 


'2n,n+l'" 


,m. 


N 


n+1 , 2n 


fin 


2n,2n 


(7.7) 


Then,  from  Equation  (2.11),  the  (i,j)th  element  of  (h^.)  is 


''ir°2''l^jii 


(7.8) 


where  jP-i|  is  the  c of act or  of  the  element  Vaai   which  is 
obtained  by  deleting  the  jth  row  and  ith  column  in  |p(  and 
multiplying  it  by  the  sign  factor  (-I)"'"  '' . 

But  this  cofactor  could  be  also  obtained  by  suppressing 

_i 
n+1  rows  and  n+1  columns  from  the  2nx2n  matrix  M   .  That  is, 

_1 
we  suppress  from  M   the  first  m  rows  and  m  columns,  the 

(n+j)th  row  and  (n+i)th  column,  and  multiply  the  sign  factor 

(-1)  '^ .   By  the  Jacobi   theorem,  however,  we  can  evaluate 

this  cofactor  by  evaluating  the  corresponding  minor  in  M  . 

From  Equation  (7.2)  we  note  that  the  corresponding  minor 

consists  of  the  first  m  rows  and  m  columns,  the  (n+i)th  row 
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and   (n+j)th  column  in  M. 
Thus 


P.il=|M 


-1 


^11' ^In'     ^Ij' 


^nl ^nn'     ^n j  ' 

-^il' -^in^-V 


(7.9) 


and,    substituting  the  above  equation  into   (7.8),   we  finally 
obtain  as  in  Equation   (2.17) 


hirVVoil 


="1-^ 


R^,....,   R^^,      R^^\ 


'  V  p  P       ' 

■    nl-""' •  *-""nn'       nj 
»  » 

-^il    '••••-^in   '"^ij" 


(7ylo) 


APPENDIX  B 
EVALUATION  OF  A  DETERMINANT 


Since  the  OS-8  operating  system  for  the  PDP-8/e  computer 
did  not  have  in  its  library  a  program  to  compute  a  determi- 
nant of  a  matrix,  a  subroutine  is  provided  to  invert  a  matrix 
or  evaluate  a  determinant  of  a  matrix.  The  Gauss -J or dan 
pivotal  elimination  method  is  employed  here. 

In  the  evaluation  of  a  determinant  by  the  pivotal  elimi- 
nation method,  two  rules  of  transformation  for  determinants 
are  involved: 

.  1 ) .  If  a  matrix  G  equals  the  produet- of  matrices  A  and  ;- 
B,  then  1C|=|AJ|B|.         .  (7.11) 

2).  If  a  matrix  B  is  formed"  by  interchanging  row  i  and 
row  k  of  A,  then  |B| =-jA| .  (7.12) 

If  a  non-singular  matrix  A  can  be  reduced  to  the 
identity  matrix  I  by  premultiplying  A  by  a  matrix  B 

BA=I,  (7.13) 

then,  by  postmultiplying  both  sides  by  A"  , 

BAA~-'-=IA"^,,  , 
or 

B=A"^,  (7.1^) 

that  is,  we  see  that  B  is  the  inverse  matrix  of  A. 

The  operation  in  Equation  (7.1.3)  could  be  carried  out 
in  n  transformation  steps 
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8^- 


SnVl---^2V 


(7.15) 


successively  diagcnalyzing  one  column  at  each  step.  For 
example,  at  the  first  stage  the  resultant  matrix  is 


E^A= 


/ 


l/a^l,0, 


,0 


-a23^/a^^,l,0, ,0 

-a^^/a^^, 0,1,0,. ..,0 


"^nl/^ll'°' 


,1 


^ll'^12"'*'^ln 


N 


^nl,^n2 ^nn 


1. 


/a^^,...,       ^in/^11 


"■iz'^-ii 


°'^22"^21^12/^11' • • ' '^2n"^21^1n/^ll 


°  •  ^n2-^nl^l2/^ll ' '  '  ' '  ^nrr^nl^ir All 


.^(1) 


(a^^?^0) 


(7.16) 


In  general,  after  k  steps,  we  want 


Ej^Ej^_^...E^A 


1, 


k   N 


•°'^l,k+l'""^l,n 


(7.17) 


k         k 
°'^'° °'^2,k+l'"V'^2,n 


'  k         k 

°»°'° ^'^k,k+l"-"^k,n 


,°'°'°'"-'^'^n,k+l"-"^n,nJ 


where  the  superscript  k  indicates  that  the  elements  belong 
to  A^   .  Then,  after  n  steps 


=A 


(k) 


EnVl-"^l^=^' 


(7.18) 
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(]r) 

In  order  to  produce  A     '$   the  matrix  E^   should  have  the 


form 


V 


1,0,0 -^ik  /^kk  •••"° 

o.i,o,...,-a^;;Vaj^;^^..  .,0 


6,0,0,...,   iAkk^.-.'"0 


^.0.0 -nk'/a|^k' i 


(-kk'^0) 


(7.19) 


We  note 


l\l=V-kk'- 


(7.20) 


The  operation  in  Equation  (7.18)  could  be  expressed  with 

recursion  formulas  which  calculate  the  elements  of  A^^ 

the  Clements  of >A^^~-^h  ,  1  - -- ■•--.--- -  ; :^ 


k   k-1/  k-1   (^k-l,  ^ 
^kj  kj  '^  kk  ,   ^^^^  ^-  ^ 


'kk 


k   k-1  „k-l  k    /■i^v\ 
^ij^^ij  -^ik^kj'   (^^^^) 


(j=l,...,n)  , 


(7.21) 


for  k=l, . . . ,n,  provided  we  define  a. .=a. .. 

The  procedure  to  obtain  the  inverse  matrix  is  then 


Let  us  denote 
,(k) 


B' 


Vk-i---V- 


(k) 


then  the  elements  of  B    are  calculated  from  the  elements 
(k-1) 


of  B' 


by  the  recursion  formulas 


(j=l,...,n)  (7.23) 
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Since,    from  Equation   (7.I8), 

i=li|=P„Vi---EiA|' 

using  the  rule  in  Equation  (7.11)  and  the  result  in  Equation 
(7.20),  v/e  obtain 


^^nn  ^n-l,n-l-- -^33  22^11^   '^l' 


or 


|A|=a^^a^2...aJJ;l.  (7.2if> 


When  the  computations  are  carried  out  by  a  digital, 
computer,  the  roundoff  error  occurs  in  each  stage  of  the 
elimination  process  and  propagates  to  the  next  stage.  An 
analysis  of  the  error  expressions  corresponding  to  the 
recursion  formulas  for  the  elimination  process  reveals  that 
this  roundoff  error  can  be  minimized  by  proper  choice  of  the 
ratio  ^j_]/\-^   (McCalla).  To  reduce  the  error,  the  numerically 
largest  element  is  selected  as  the  element  a,  ,  ,  which  is 
called  a  pivot  element,  from  the  elements  a.  .  (i,j  _^k).  In 

J 

order  to  bring  the  largest  element  to  the  pivotal  position, 
rov.'  and  column  exchanges  are  necessary  and  the  corresponding 
sign  change  in  the  determinant  should  be  applied  by  the  rule 
in  Equation  (7.12). 


APPENDIX  C 
NUMERICAL  INTEGRATIONS 


1).  The  error  function  in  Equation  (3.13) 

1  X        2.  (7.25) 

erf(x)=2Tf"^  XcLtexp(-t  ) 

0 
i.  evaluated  by  the  rational  approximation  (Abramowit.) 

errU)=l-(a,t.a,t^.a3t3.a,t  V'-'^'-''^'^' •  '^•''' 

(0^<°«) 
where 

t=(l+px)"'^s  ,  .  ,  r 

p=o. 3275911.  ^ 

a^=0. 25^8296, 
a2=-0.28i^'^967. 
a3=l. ^21^137, 
ai^=-l.  ^53152, 
a^=l. 061^05. 

The  error  bound  in  this  expansion  is 
|r(x)l^  1.5x10"'^. 
2).  To  evaluate  the  statistical  function 

-i  1:  ,    i4-2x  (7.27) 

Q(x)  =  (2tt)  ^;dtexp(-2t  ), 

X 

the  following  expansion  is  used  (Abramowitz), 
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.    ,  88 

Q(x)=Z(x)(t)^t+b2t^+b^t^+b^tVb^t^)+r(x),        (7.28) 
where  (0^<=o) 

Z(x)=(2-rT)"*exp(-ix^), 
t=(l+px)~^, 

p=o. 2316^19, 

b^=0.31938l5, 
b2=-0. 3565638, 
b^=l. 7814779, 
b^=-l. 821256, 
b^=l. 330274. 

The  error  bound  is 

|r(x)|^  7.5x10"^. 

For  a  large  value  of  x,  the  approximation  of  Q(x)  by 
continued  fractions 

Q(x)=Z(x)[l/(x+l/(x+2/(x+3/(x+. ..)))],   (x>0)    (7.29) 
would  be  a  better  evaluation.  But,  since  the  difference 
between  the  two  approximations  was  not  significant  enough  to  show 
up  on  the  plotted  graphs,  the  expansion  in  Equation  (7.28)  is 
used  in  the  computation  throughout  the  range  of  x. 

3).  For  other  numerical  integrations,  the  trapezoidal 
rule  is  employed. 

Suppose  f (x)  is  defined  at  equally  spaced  points  x. 

(i=0,l, . . . ,n),  then  the  area  under  the  curve  of  f(x)  between 

Xq  and  x^  is  computed  by  the  trapezoidal  rule 
n-1 
%y^^^K)-^^(^i+h)]'  (7.30) 
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where  h  is  the  distance  between  two  points  of  x 

^"^i+l'^i*   (i=0,l....fn-l).  (7.31) 

Writing  the  above  formula  in  expanded  form  and  substi- 
tuting the  relation  x.  >=x.+h,  we  obtain  the  classical  form 
of  the  trapezoidal  rule 

.   T=-|h[f(xQ)+2f(x^)+2f(x2)  +  ...+2f(x^_^)+f(x^)].     (7.32) 


APPENDIX  D 
INTERPOLATION  BY  SPLINE  FUNCTION 


When  an  interpolation  of  data  is  needed,  the  spline 
function  is  used  to  approximate  the  intermediate  values  of 
the  data.  The  spline  function  is  a  numerical  representation 
of  the  curve  of  a  flexible  strip  or  spline  used  in  drafting 
to  draw  a  smooth  line  through  a  set  of  points  (Carnahan) . 

Suppose  we  have  n  functional  values  (x.,f(x.)) 
(i=l,...,n),  and  want  to  interpolate. the  data  over  the 
interval  (x^,x  ).  The  restraints  that  define  the  spline 
function  over  interval  (x^  »x„)  are  then  as  follows: 
1).  In  each  interval  the  spline  function  is  a  cubic  poly- 
nomial, that  is 

p^^^(x),  (i=l,2,...,n-l)  (7.33) 

where  the  subscript  in  p„  indicate  that  p  is  a  polynomial  of 

degree  3- 

2).  Each  cubic  polynomial  p„  . (x)  should  match  the  functional 

values  at  its  tv/o  ends  x.  and  x.  ,.: 

1      1+1 


P3^i(x.)=f(x.), 


(1=1,2 n-1)  (?.3^) 


P3,iK-M)=^K4-lV'  J 
3).  The  first  and  second  derivatives  of  successive  polynomials 
must  match  each  other  in  value  at  the  intermediate  points: 
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P3.iW)=P3.i-lK)' 
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(1=2,3,. ..,n-l),        (7.35) 


and  the  second  derivatives  vanish  at  two  end  points: 

P5,l(x,)=0. 

(7.36) 

p"    , (x  )=0. 
.  •^3,n-l  n' 

Equations  (7.3^)  through  (7.36)  amount  to  4(n-l) 
relations,  which  enables  us  to  determine  the  ^(n-1)  coeffi- 
cients needed  to  completely  describe  the  spline  function  in 
Equation  (7.33). 

For  convenience,  define 

^i"^i+l~^i'  (7.37) 

ei=P"3,i-i(^i)=P"3,i(^i)- 

Since  each  p"  ■ (x)  is  a  linear  function  of  x,  it  can  be 
obtained  on  the  interval  (x- ,x._,_.  )  by  linear  interpolation 
of  the  values  g.  and  g. ,.  at  each  end  of  the  interval: 

P3,i(x)=g.(x.^^-x)/h.+g.^^(x-x.)/h..  (7.38) 

By  integrating  twice  and  imposing  the  conditions  in  Equation 
(7.3^),   we  find  that 

P3,i^''^=Si^'^i+l-^^'^/(^N)+gi+l(^-^i)^/(^^H^  + 

(:fi^lAi-gi^lV6)(x-x.)+  (7.39) 
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if^/h^-g^h^/6){x^^^-x),      (i=l,2,...,n-l). 

By  differentiating  Equation  (709)  and  imposing  the 
conditions  in  Equation  (7.35)i  we  obtain  the  following  system 
of  linear  equations  in  the  g. : 

gi.lhi_iAi-<-2g.(l-i-h._^Ai)+gi+l 

=6C(f.^,-f,)/h,-(f,-f,_,)/h,_,]/h,.       (7.^0) 


(i=2,3,...,n-l), 


and 


gl=0, 


gn=0,  (7.^1) 

which  follow  from  Equation  (7.36). 

Note  that  Equations  (7.^0)  and  (7.^1 )' amount  to  n 
simultaneous  linear  equations  in  the  n  unknowns  g.,..,,g. 
Thus  the  values  of  g^ » . . . ,g  can.  be  determined,  and  the 
substitution  of  the  obtained  g.  (i=l,2, . . . ,n)  into  Equation 
(7.39)  yields  the  individual  cubic  spline  polynomials  for 
use  over  successive  intervals. 


APPENDIX  E 
LIST  OF  PROGRAjyiS 


9^ 

PROGRflH   L0H6G 


COKKOU    P,ft/U,U/nLP/ES.MUG.T,PlvUP.SIGKn 

DIKTUSION    P<6,6),flC6/6)..LICe.,t.).U<6.6),ftLP<3>.ESC3),l';UG<:JGe) 

PIKTHSICM    TRLK200),DftTn<20e> 

P1^^3.  I41592£.536 

tM:!.Tt:Ci.  1 

PD    10    1=1,100 

10  kl:o':i>  =  -13?& 

DO    100    HRXT^i, ISO 
T^r-n.T'^Fl.ORKMnXT? 

[)UKO  =  R-C 

P:D=="t)U{{0-'T 

Rr'l3==2.0+C>LiUG''T**2-R  -  * 

Re=-i.e 
Rr)e=:0.o 

RDDO-- 1.0/3.0 

P<J ,  I>  =  RO 
P<I  .2>  =  R 
f-'<Ij3>-Ri'6 
P<l,4>  =  Rr> 
P<2,2>=R0 
P<2,3>=-RD 
P-:2,H?  =  Ri>0 
F<3,3>  =  -RD[)0 
P<3/4>  =  -RD[) 
p<4,4)=.~R[)t)0 

DO    200    1=2,4 

DO    260    .1=1,  MUD 
200  P<3  ,  J)  =  P':J,  I  > 

SUDl^l.e-R*';2 

DO    250    1"1>2 
DO    250    J  =-1  ,2 

DO    235    ll=--l,'l 

DO    235    Jl=l,4 
235  ft-:)  I,.U>  =  P<1  1  ,  Jl> 

DO    240    L^-I,4 
240  R':3,L)  =  R<:2+1 /L) 

DO    245    L=l,4 
245  fi<l  /3>=fl<L,2+J> 

KRT=^^3 

CRI.  L    lira  I  »    <fl,llRT/D> 
250  U<3  ,  J)  =  D''SUDI 
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PO    2701    1  =  1,2 
[)0    27D    J=I,2 
2?Q  y<)vJ)  =  U<I,J)y'<U<l,  I>*U<J,  J>)**0.5 

t>UKC=<U<I  ,  1  >'«:U<2,2)>**0,  f.;/'4.0-'PI*  +  2-^<1.0-R**2)*>J(>.5 
DUM-d  .  e- UCI  ,2)**2>*i:0.  ti/Us  j /2> 

H  <f)UI'. .1)275,  2£:0,  250 
27r.  Rrvt-F-l+flRC 

2&6  Wpr:  =  r)UMC'«<  <  1  .  R-y<  1 ,  2)*r2>*  +  0:  5-y<  J  ,  2>*fiRC> 

yP"SQr:T<-RDDa>/'2. 0^pi 


P<J,I>=R0 

P<3,3)--R 

f'<i  ,4>  =  RDG 

P<},6>  =  Rt> 

P<2,2>=Ra 

P<2,5>  =  Rt)e  - 

P<?,3)=R0 

P<?,'1)  =  -R[> 

P<3,fe>  =  Rr)0 

P<'1/^)  =  -RDD0 

P<H,6)=-RDt> 

P<t;/5)  =  --RDP0 

P<6,6)=-RDD0  ' 

Cftl.L    L0K6I 

PCf7  =  WPri/UP-SI6MR 

TRL;<MfiXT>  =  T 
le.O  [)RTft<:MRXT>  =  POT 

CftUL    OOPEH'C'SVS'^'ZXiiO') 

l.MMTE<4,?0O>    HfiXT,  CTRLK]  >,DRTR<I>,  1  =  1  ,IiaXT> 
700  FDr<{5RT<fti',40OR6> 

CRI-L    OCLiiSE 

EHD 
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f.Ul.kOUTJME  LOUGl 


[JJJU.HSICIN    p<6.e>/fi<6,6).U<6.6),  U<6/6)/RLP<3)VES<3>/l1UG<JC<e) 
DS)C>-C.O 

G-J.O 


5PP 

rj'3C'' 


Tf-G^DPELT 

1  F  <  T-TF-DTF  >2fi0,  200,  21  0 

£:3Gi;R:^f.:iGl'.rwl)i:lG 

GO    TO    ifJO 

1  F  <  G-  2 .  e.  1  >  226 ,228 ,  23D 

S:iGi;n=^S)GI'.R+I>SIG 

GO    7  0    36!:. 

«:)  GI'.R==S1  GMR  +  2.  6*DSI  G 

TS:'-T-TF 
S==i:l!KTF> 
C:'=(;as<TF> 
FiF-=S/TF 

RFI)---DUM0/'7F 

r:r!>!>^2.  t"'Dune:''TF*^^2-RF 

C"CO&<TS) 


P<5 ,2>=RF 
P<}  ,S>==RFD 
P<2,3>=^RS 
P<2,'i)  =  -RFt> 
P<2/«^.>=-RSD 
P\'3,5)  =  -R£E> 
P<'(,5>  =  -RFt)t) 
P<5.,d.>  =  -RSDt> 


3J0 


DO  310  1^2,6 

DC  35  0  J=l  ,\\U[) 
P<]  ,  J>'=P<.T,  1> 


320 


DO  320  1-1,3 
DO  320  J =1,3 
Ff:)  ,  J>  =  P<I,  J) 


!inT^:3 

CftLL    l5nTlM<n,ltnT,D> 

SUl>2^-D 
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DO  3li0  1  =  1,3 
[>Ci  350  si  =1,3 
DO    33^;    n  =  l,G 

DO  33r:.  Ji  =  i  ,e. 

33?  f{<)  l,Jl>  =  F<n,  Jl> 

DO    3<ie    L=l,6 
3<10  n':'i,L)  =  fK3+I,L) 

DO    3^5    L=:l,6 

3'U.         na.  ,''.>=R<L,3+j> 

CfiLL    MRTI»<R,nRT,D) 
3'c.t^  LK)  ,  Jy-D/SIJD2 


3?0 


DO  3?0  1=1,3 
DO  3?0  J=I ,3 
kK  )  ,  J  >  =  U<  1  ,  J  )  v<U<  1  ,  I  >*LKa  ,  J  )  >*'J  6.  5 


DO  3&?  I0=i,3 

MUl:=--ICt+3 

ri'^lF-'EiKKUD-'S)*! 

J2=^ir;EM<MUD-'3>  +  i 


Dun 

DIM; 

1F< 

301 

DU11 

Dur. 

ir< 

302 

Kil.P 

Dun 

DUI' 

Dur. 

]f  < 

30'{ 

DUIt 

f{R(.: 

1F< 

3oe< 

RRt 

SOl'i 

ES-: 

60 

50? 

ru- 

38f. 

es:-: 

307 

coil 

0=-U<I2,  I0> 

6-nBS.<:i>uHG>,- -    ■■  ;^-^- -.   ^ 

DUMO-Ci.  9?(>381  ,333,  333 

c^=uau,in 
ejr^ftr::s<DUi;o> 

DUUO-e.  99>382,  303:.  3S3 

<I0)  =  U<I2,  IO>+U<Kt..  Ii>  +  U<n  ,  I2> 

C:^U<12,  I0>*U<10,  n>-U(II,  12) 

0r.[)UM0/<<l.C-U<I2,  I0>H:*2>*<1.  e-U<I0,  II>**2>>'!»J0,5 

l^RE;i:a>UI5Cf> 

DUH 1 -■  I  .  0 >304  ,  3.86 ,  306 

e:--<  }  .  G-DUIUJ*4=2)*H:0.  S^'DUtlD 

=-^flTRK<:DL5|iCO 

DUiiCO 300  .305,385 

=  Pl  +  RFiC 

ICO  =  RFcC 

TO    38? 

(10>=0.0 

!0)=^0.0 

TIliUE 


|;R7  =  3 

CRl  I.    MRTIlF<V,nRT,D) 

D=-P.OO<D> 


DIUU^--<U<1  ,  i  >i-ll<:2>2>*LI<3,3))-«-*8.  ! 
DUIH-D*:-*©.  livE5<l  )  +  RI.P<  J  )KES<2)- 


PJ  >  +ftLPC2>  +  <ES<3)-PI  )^>ni:p<3,v 
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t>S3  0=:iJprjtt-^up>;r)[)ELT/'2.e 

G"(vH.O 
GD    TO    3e0 


100 


RfifURU 
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THIS    PROGRRfl    IHUCRTS    R    IIRTRIX    Rv'H.IO.     THE    OETERniK'RHT    t> 
or    THf{   isRTRix  r<:n.n>    is  RLSO  CRLCULRTED. 


l>IKf;HSIOH    R<6,6),  IP<6).  JP<:6> 


IP 


?5 


2!i 
-I? 


5P 


('l:> 


?tv 


m 


DO  ie  i=K,n 

DO    20    J=^K,H 

DUI!1-:R<I,  J> 

DUP.l^RBSkDUHn 

1  r  <DUfU  -DLinO>26,  28..  25 

IP'CK):^! 

JP<K>=J 
DUI'0^r[)LIMl       ' 

coiniwuf 

JF<DUIia>S!20^9Je,45 
1  r-  <  Jr- < K>-K)5i20 .  52 .  '1? 

DO  r;o   1  =  1  ,H 

DUIUl--R<l..K> 

R<)  .K>=^R<1.  JP<IO> 

R<) . JP<K>>=DUMQ 

1 F  < 1 P  < K > -K > 926 , 65 , 54 

DO    55    .1=1  ,H 

DUnC.=^-fi<K.  J> 

R<K. J>=R<IP<K). J) 

R'CJF-viO..  J>  =  DllHO 

PlU-R<i:,K> 

D'^DiplU 

DO    ?0    I  ==  i  .  H 

1F<  J-:K)75,?0.?5 

R<)  .K)  =  -R<I,1C>/'PIU 

DO    ?0    J=1,H 

jF<j-K>eo,?o,eo 

R<  )  .  J>  =  Rn  ^  J  >i  R<  1  ,  IO*R<K,  J  ) 

cojrniRiE 
DO  r^e  j=^i,M 

R-CK,  J>  =  fl<i<,  .Oz-plU 
R<K..K>  =  1  .O/'PIU 

H  <K-IOI0,  10,  I  15 

HM'-N-l 

DO    130   f!=^l  ,iin 

IF<JP<I>-l>92ei,  127,125 


100 


J  25 


126 
12? 
1 3V. 


136 
9  IP 

«^f2t« 
?21 


DO 

120    J 

=  1,H 

DUi'o=--n<J 

,  1> 

ft<J 

.  !>  =  -- 

H<J,IP<i)> 

R<J 

. 1P<I)>=DUM0 

COP. 

1 1  {{i.ir: 

ir  < 

Jf'CJ  > 

-1>920. 138 

■  135 

DO 

130    J 

=  !.!< 

Dur 

tir.-fKl 

.  J> 

ftO 

...1>=- 

R<jr<)  ),j) 

R<..U-<1>, 

.n  =  DllM6 

COl,' 

riiiijE 

60 

10    99 

9 

l^Ri 

T[:<i. 

9I1> 

FOM'ifrK' 

HRTRIX    IS 

SilJGULRR.  '> 

GO 

10    99 

9 

I'K  I 

TF.a, 

921  > 

FOK 

?tfn<' 

HnROURRE    ERROR. '> 

fiHi 

UR!{ 

ENI) 
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PKCiGRfiM    STOSK 


PlfifiKSICiM    •rRiJ<20O).[)ftrft<26C>,SH<5O) 
KR}sl--5 

inc}=i 

r.RK2-=9 
1HC2-1 
S:C:f;[I=-I.D 

sent  2=e. i 
[)R==e.o 

DO    100    LF'i:^inin,HRJ<I.  I»Ci 
R"i:CRLl^FLORT<LPl> 

DO    iOO    LP2=-HlM2y  HaH2/lMC2 
C=-i:CRL2>:FL0fiT<LP2> 

CRl-L    SI  HO 


CLR 

CLL 

TRP 

•vLPl 

TRfJ 

<60 

RTI. 

RTLiRTL 

TRD 

<5& 

Dcn 

NDR 

CLL 

TRD 

NLP2 

TRr> 

<e.e 

R7L 

RTL^RTL 

DCR 

VDR<* 

CffLL    OOPEIK'SVS'.PR) 

i'Rnr<'!i,2eio>  nDR..R..Q,  <TRua  >>drtr<i>>i==^i.kdr> 

260  F0kHRT<ft2/'ie.2fi6) 

CRLL    OCLOSE 

100  coirriMUE 


102 


SUI^kOLITIME    S1H0 


COnnOU    HDR,R..a.TRU,DRTR 

P 1  HEWS  I  ON    7  nU  <  200 ) / DRTR  < 200 ) 


0.  25'1S25»59 

-0.284496?'! 

1.42MI37 

-  I.453IV.20 

$.0614  604 

0. 327591 1 

1=26 

3. 1415926 

=  F-l''FLCifST<NDr-n 

^Pl/'4.6 

(V2.G:/PI 
2-^2*WDPl  +  J 
i;i.f  =  3"H[)Pl.'2+i 


El=^ 

V.'S- 
F.4- 

E6=^ 
P3  5 

Pi: 


DO    200    LP I -5, 150/5 

T--C.  l^FLORKLPl) 

T2»a';:H:2 

t>UliCi:-O«Tv2.0 

C:a=^C:DS<DUMO> 

S(^=--SJiKDUHO> 

R'-?11KT).''T 
l>UIUt=R-COS<T> 
RlJ=--DLinO/'T 
Rr'I:^=2.e>![)lJHC''T2-R 

RrM3  0=-i.o-'-3.e 


C:0=^^KD 

C5-I.0-R 

C:2'-I.0+R 

C:3=--  RDDO+RDD 

C:4=^^^-RD[)0-R[>D 

C5=  C1>;C4-C0*':2 

t:6--C2-^C3--C:0*>:2 

C:  V  ^-  C.  0  *  C  a  +  0  >^  C  2  V  S  Gi 

CV-  CO-<f.O+  Q:*C  i  >:-C0 

[)I:T=-C5*C6 

>:o=-ci^:C6 

Xi-^C2«C5 

X2-  <R*C7 >*r=.2^C2/'C6 

X3==<fl^C!E!)'«:*2/(;l/'C5 

X4  "  <  R«  C0.''C2  >  =<=* 2 

Xf.-  ^R'^-Sa^'Ci  >**2 

X6-X2+X4 

XV^>;3+X5 

>:8"<(;5*C.6>-«K.0.5 
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[)0    100    LF-2=I/nLP2 

DUMO- FLCiriT<!'.iJt)e> 
TI!i:=^DLI!50Hpn-pi 
CT'=COSaHE> 
ST'  SIIKTIIE? 


[)3:=X6«C:T*  +  2+X7*ST>i=+-2 

Xfi=^Cr.*C?+CT 

y.iV'-cf.'^c.-^iSi 

lF<D3-?3.e>46..47.47 

-If. 

EK{)3--=EXPC-[)3> 

GO    TO    1f>' 

17' 

EX!>3-0.0 

19 

sun 2- 6.0 

PO    IGO    LP3==1  >nLP3 


l',ff:C=LP3-l 

l)UIUV=FLOriT<KLiDO> 

PHI==DLI.Ma*Pi  3 -PI  2 

CP=^C:CiS<PH!) 

SP=-S3!KPHI) 

t)lll'0=-2.e'iPHI 

c:2P=cos<:t)iino> 


t)i:=XI>?CP»i*2+X0*SP**2 

l>2---fi*<X9>rC:p+XiO*'o:p>/>:8 

D1^-D2/DiK:;6.5 

1>8-D1'i>:2 

[)l>=-D3-[)8  . 

E0=^  i  .  0.''<  1  .  Ct  +  Et'.*DLiril  > 

[)l!riC^-El<E6  +  E2:=EO**2+E3*E0*>i'3  +  E4  +  E6»i*4  +  E5*E0**i; 

)F<D8-?3.  e>t>f?..  51  .tu 
DO  EPF=  1  .  0-[)UUa*EXP<-DS) 

lF<D1>L.'1,55,tV5 
li4  EPF^^-EF^F 

GO   TO   e.5 
SI  ERF^l.O 


i:.t"i  Dl!P.O=<!  .  C'.  +  [>5>H;EX[)3 

60  3F<r)t.-?3.  e.)<-..I  ..  62ji62 

61  ijnfn^^pi3*[)-ie<:i.5^[)S>i=ci.  o-erf)*expc-d5) 
GO   TO  e^j 

62  [>uin-o.o 

Co  sun- <Duno-PUin  )HC2P/'[>i*^^? 


10^ 


ir<LF*3'-l>70/70,72 

?c     sun2=suM2+sun 

GO    TO    150 
72  ir  <LP3- ML P3>7.'.  ,70,70 

74  j:ui:2-sui'.:''+2.  e.^?.mi 

H:.C!  COHTIHUE 

lF<LP2-nS0..  80,82 
80  8l.!ni"SUIJl+<;:LlH2 

GO, TO    iOO 
82  ir  <LP2-HLP?)e.4.,80,8e 

8-1  SiUliS  ^SLIIII  +2.  e=*3UJ12 

JOG  CONTINUE 

l'rM=^DET>;*I  .  r.*?-UI11*Pn«H:2-^PI<=*3/32.  0 

iri=^ypn/fO 

K!)fi  =  NDR+l  - 

TnU<HDn)=^e..  i*FL0RT<LPi> 
2C-tO         •     t)ftTF{<HDR)  =  tJT 

P.S^TMRN 
ENJ> 
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PRCiGK-ftK    COREi: 


cciMiJOM  iU)n,R.fi>TRU,[)fan 

D 1  nr; ns i  on  7 ru < 260  > ,  drt r  < 2tic  > ,  mug <  1 00  > 

r{RXl=5 

r.RX?---? 
sen  1. 1  =  1.0 

$:(:fiL2---e.  1 

c>Tnu=o.5 

CfiV'O^-  <  1 .  0.'3 .  Ci >  *^Q .5 
C:(il<  1  - 1  .  0/  <  2 .  e* P  n * *Q .  5  ' 

t>{:i.T"pi/'2e.c  : 

fTI -0.  25'1S2959 
E2-^-e.2t:44  96?.4 
E3=-i. 425413? 
E4^-i.453I52e 
EU.:- 1.06 14  05  4 
E  6^^  0.3275911 


DO    t>0    1  =  1  ,  100 
50  nLlG<l>  =  -I3?e. 

DO   100  LPi=nini,itRXi,  iiici 

DO  100  LP2--J!IM2.HfiX2.  inC.2 

S  Cl-R  CLL 

S  TRD  \LP1 

f:  TRD  CGO 

S  RTL;PTLiRTL 

f:  TRD  <:p&^ 

S  DCn  XDR 

5  Cll- 

?.        TRP  \LP2 

6  TRD  <6e 

5;        RTL  ;RTLif;TL 
S        DCR  \DR« 

CRIl-  lOPEIK'^-.VS'.DR) 

Rr;nD<4,2Liio)  iu)R,r.o^  <ctru<j  )>DnTfici>,  i  =  i /HDn> 

200  F0RIJRT<fl2^4«2P.6) 
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3e.c 


FOrdtRK'I'IOX..  ''fl='F3.e,  iOX,  'O^'F-I.  !>'''/'> 


DC    326    IM=J  ,nDf\  I 

32ei  I'i,.  n  E  <  3 ,  32 1  >    T  fiU  <  I H  >  ,  DftT  fj  <  1  H  >  .  <  HUG  <  J  >  ,  J=  1  .  HUDO  > 

321  FCHj;faAF?.2.f:?6.?,t;x,  ioop.i> 

l'RJTF<5,323>  " 

323  FOrclSilT  <//'///'> 

sujvi=o.e 

[)0    4eO    LF'5=-1,2I 

KUi:0=LF-5-l  > 

7HU-Cl>::FLCiftT<HUt)0)^26.  Ci ' 

£:=--J:Jl(<THF) 


AJfil. 

Dun 

ECt^^ 

t)Ui< 

[iur. 

1F< 

'Jt.C 

ERF 

1F< 

l^.-l 

E  HI- 

GO 

'ir.i 

ERF 

-iji? 

DIM' 

0"E:';=S/2.  6**0.  5 

l-nEJSXUP.LO> 

i.  0/<I  .  C'.•^E6*0UlU  >> ,  I 

C=^-E  l*Et^^E2*Et1'»  *2  ^E3*E0**3+E4i  E0>«  *4  +  E5*EG-i--i-5 

l==URL0  +  -;2 

DUIll-73.  C<>4S0,45i.. '15i 

=  I.C-DLIfiO*EXP<-tJUlil?- 

U  R  1. 0  )  4  5  4  •  '1  !r.  5  i  4  5  5 

=-ERF  , 

TO    455 

=^  I  .  C 


l>l!l'0^:-CR-^C>**2/'2.8 
URL B^-COH i »;  EXP < DUISG > 
[)U/1Ct--<E:';=S>**2/'2.  6 
UPi!.l---C0NivEHP<DUnGO 
[>SU!5i=UftLO*<:MnLl+B+S--'2.e*ERr) 


lF<lRei-I>4t.l  .46i.462 

GO    TO    460 

462  IF<LP5-2I  >4*:-3.  4f.i  .  461 

4 63  sun  I  =-SUK  1  +  2 .  £;>;  DSUH 1 


4  00 


coinmuE 


l.'P'-POSE>!SUISi 


=  r)ELT/'2.e 


sun2-o.o 


X>0    'J2&    1  =  1,IU»R 
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554-; 


I  r  <  !'.Uf.)C'.- «<e  ) 596 .  590 .  59 i 

r.uf;  0:^90 

CCUniNUE 


1F<1-1  >5IC,5K<.512 

51  e  sun2^sui{2+[)fnn<i> 

GO    TO    52C1 
5i2  IF <r-HDR>5i 4, 510,510 

SJ^S  Slin:-f==SIJ!t2-*2.  OtDfiTRC  l  > 


526 
521 


I'H; 3  TE  <3 ,  521  >    TFiU <  I  >  ,  DRrfi  (  J  >  .  CHUG <  J  >  ,  J  =  I  ^  HUDO  ) 
F0r;!5RT<r?.2.E26.  ?,5J<,  I00ftl> 


536 


friOE;=0,5>!SI.n{2';:DTRU 

Forjim  <:/'/'/^iOX'«.''<:2*pi  j^'f?.  4*  iox'up='f7. 4,  lox 

'PKOBfit:ILl  rV='F5.2> 


3  OP 


COK'TIKUE 
EMI) 
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PROGRRM  RETl 


t)  I  nth'S  1  OH  T  RU  (  200  >  ,  DRTR  (  260  )  ,  MUG  <  i  00  ) 
ORl.l.  IOPE:m<'S:VS', 'ZXHO 

Rr:fiD<'i,2e>  im>r.  <TRi.Ki>,i)R7n<i),  1---1  ,hdr> 

2i£.  F0fv!{RT<R2,400R6> 

[)0    38    1  =  1,  100 
SZ  l'.Lir3<l  )  =  -13?6 

DO    ^0    1=] ,MDn 

[>uno-2eo.  o*L>nTfl<i ) 

rrjl:Or=]F  lX<DUIi0> 
10  URJTECS.-ll  >    TRlI<I>,r)R7RCl  >,  <:MUG<J>,  J=1,I1UD0) 

1i  F0rd{RT<F10.2,E;26.  ?,5X,  iOO<fli)> 

El«> 
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PKOGRRH  RET2 


COHMCm    Nr,)R,n/(v'/Tftl.l»DftTn 

Dll'.tHSIOU    TRLK20tO,DflTR<7£iC>,MUG<iOCi> 


r.iin  =  i 

IJRKl^S 

IMC  1^=1 

n3M2^^3 

iinj<2:-9 

JI<C2=-1 

S(fiLj  =  l 

0 

s(.nL2^e 

1 

DR==e.O 

t)0    1.0    1  = 

=  1, 

lec 

i'ur:a>^- 

-13 

76 

DO  ICO  LPj^rnin.MflXiviMci 

t)0    iOO    LP2=^inM2,MRX2.  1HC2 


CLR 

CLL 

TRD 

\Lfl 

TRP 

<60 

RTl. 

RTLiftTL 

TRD 

<^.6 

[>c-n 

\DR 

CLL 

TRP 

\LP2 

TRP 

<eo 

RTI. 

RTLiF:TL 

PCR 

\DR* 

CRLL    lOPENCSVS'.DR) 

fi\-:UlHA,7.60>    l«)fi,R,C!,  <TRLKI).DfiTR<l  >,  I  =  I,HDR> 
200  F0RI5RT  <R2,4e2RG> 

lM5]TL<3,3eO>    R>D 
300  FOIUlRT<'i''IOX/  'R^'FS.O,  lOX,  'O.^'F-I.  I''/'/? 

SUK^O.O 

DO    326    1M=1.WDR 

PUI!0=:i'06 :  Ci*f)RT R <  I M > 

laiDGr.lF  !X<[)UrtO) 

]  F  <  IJUPO-  90  >  30^j  >  305 .  306 

306  nuivc--90 

30J^  COinJM'JE 

sun- SUIMl)RTR<?.H> 

320  I'R  1  TF  <3  ,  321  >    TRU  C  I  M>  ,  DRTR  C  I  K>  ,  CHUG  <  J  >  ,  J!=  1  ,  liUOO  > 

321  FOIvltRKF  ?.  2.F20.  ?,5X,  IOO<Ri  >> 

1'P1TL"<3,35:.0>    PROB 
350  F DMtRT  <  /'/'-^  I  OX ..  ' PROBRB I  L  1  T  V^  "  F5 .  3  > 

ICO  ccurriMUF 

EN!) 
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fKOGRftK    [)P0L 


OIKKMSIOIJ    TRU<36),f)RTH<3ei>,MUG<]0CO,DRP0(31>.TnP0':3i>,PRr><31> 


DO  r.o  1=^1  / 100 
tie  ni.io<i  )=-i37G 

[>C    iCO    LPl-3.5 
DO    ICC   LP2'-3,9 


s 

Cl..n    CLL 

i: 

Tfir:-    \LP1 

f: 

Tf'.P  <eo 

S 

RTI  .jRTLjftTL 

s 

TRr:     <:5<:. 

t; 

DCf{    SDR 

S 

CLL 

f; 

Tftr:-    \LP2 

f: 

TP.r.     <63 

S 

RfL;RTL;PTL 

s 

DCn    \OBM 

199 

tfJLl.    IOPEK< 

RLf;tJ<. 1,260) 

a-tto 

F0MinT<ft2,C.: 

IU'5R,R.Q,  <TRLKI  >,DRTfK]  ).  I  =  1,NDR) 


VR]iTK<3,3eO>    HiO 
3PC!  F0Rlin7<'l'I0>s, 'n='F3.6..  lOX, '0='F4.  l''v/> 

Se-i  l*R]Tn<l  ,3C-i5>    R,Q 

3or.  Forv!inr<2nc.  i,  iON'MRxiriUii  TfiLi=='> 

Rr;Rt)<i,3iC>    TU 
310  FORf5RlCF20.  2> 

DLIIiO^^rrt/O.  J:-.-vO.  61 

np=- iF]x<Duno> 

iHi5  2-^e 
i»i>i^np-2 

DO    -'.OO    LP3:-IHt)l,NP 
I){!>2"IM()2+1 

Tnr  o<nu)2>=:Tnu<LP3) 

•lOO  DRP0<IH[)2>=:DRTR>CLP3> 

■rftpo<'i>=r5.  0 

DRP0<'J>  =  6.0 

H3=^30.C-FL0RT<HP> 

02" 1.0 

C3=-  1 .  C-'HS 

C'i- 2.  0*<1  .  0+1  .  Ox'HS) 

Fi'=^C.0t<DRP0C3>-DRP0<2>-DJ"iP0<:2)  +  [)RP0Ci  )> 

F2'^ 6 .  0^'H3*  <  <DnPO  < 4  > -DRPO  <  5 >  )  -'H3-DRPCK3  )  +DRPO <2>  > 

Pi"C5.0 

P2-=<Fi>;:C4-F2)/'<Cl*C'1-C3> 

P3=!i"Cl«P2 

P'J^-fi.O 


Ill 


699 

HOL 

00 

r-ci^: 

[)UI'. 

[>un 

l>U."i 

jr  < 

i.2& 

ir  < 

6Ki 

Pc.=- 

00 

C.3P 

HDL 

cm^ 

[)RT 

t>-^  100.0 

600    LF"s  =  NP.38 

=^FLOnT<LFM) 

I~F'3/C..  0/n3*<36.  0-P0S>*>='3 

;:=^<[)RF'0<3>''H3-H3^P3/f.,  0)*<:38.  O-POS) 

o-r>uin+i:)ij!'.2 

l)LI!lO>eiO,i;.28,628 

f!OL[)-[)UiIO  :■  e  1  0 ,  630 ,  638 

0.9i.'.'^P3 

TO    6<^'3 

fJ^OLIMO 

n<:i.p-4)=Duno 


siwt^^o.o 

DO    tiOC    LP5==I..  38 


J/06 

t".0^ 

boo 
ivio 


SI  5 


DUnO-^100.0*DRTfKLP5> 

ivJi:o=-inxa>uiio> 

ir  <MLiD0-S0)50t">.5C5/506 

l',UJ;0=9G 

$:UP;=^SLI(l4[)RTn<LP5) 

l'Rr!E:<3,510>    TnU<LP5>..DRTfKl.P5>,  <Mlfrj<i>,  I  =  I,J-}UDO) 

f0kl'iRT<r-7.2,F:2e.7,5X^  iOOfiO 


PR0r:-SU}5vG.  Z>    ■ 

l"R]TE<3/5I5)    PROB 

FOMSRT  < yyy  I  ON  ' PROBRB I  L I T V-  "F? .  4//'x'/> 


1'P1TF<1  ,700> 
?00  FOMJRTC'RPPROUE?    THEN    TVPEHl    !.'> 

FUIfUKi  ,705)     in 
70r.  FOIvKRKin 

!F<lH-i  > I  99, 728 > 199 
720  CBI.l.    OOPENCSVS'^Oft) 

UR3Tn<4,730>    NOR.  fl.  Q,  <TRIJ  <  I  >  ,  DRTRC  I  )  ,  I -1 ,  NDR) 
730  F0RKRT':R2,62H6> 

CRI-L    OCLOSE 


100 


CONTINUE 
ENP 
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PROGRftM  SfiPOL 


ccip;rjon  tfilic<.. [>infsi:i.TftLi,[)ftTn 

t> 3  J'.f: HS  1  OH    1  RtJO  <  1  50  >  ,  DHT f{t3  <  I  50  >  ,  MUG  <  1  CO  >  >  TftU  < «  ,  30  ) 
OintHSIOW    DRTfKS,  30>  ..HPCK6)  ,  [)HP0<6>  ,  t.0P0<6>  ,  PHD  <6.> 


50 


KP0=^3 

DP    '-.O    I=-l,  100 
nU(^<l)  =  -13?6 


JOO 


CHI  I.    I  OPEN  C  '  SVS  '  ,  ' ZXi,'0  '  > 

R(:fiD<'1 ,  1  00  >    lUJR  ,  «:TfiLiO<  I  >  ,  ORTROC  I  )  ,  I  -  I  .  KOR> 

FOI<I5RT<R2.^0OR6) 
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1>D    J  20    1  =  1,30 

TRIKI  ,  1  >  =  TRUOCMUC>0> 
t>RTR<  1  ,  I  >  =  f)Rffl6<IJLI[)G) 


DO    200    LPi-1,9 

DO    300    LP2-l>r. 

1HIJ  =  LP2+I 

S  CLR    CLL 

S  TRD    NLP2 

f:  TfiP    <60 

S  Fm.  JRTL^RTL 

J:  TRf>    (S*'. 

S  DCR    \Dfi 

S  CLL 

ft  TRP    \LPi 

ft  TRO    <60 

ft  RTLiPTL;PTL 

ft  DCR    \DRis 

CRLL     10PEN<':ftVS'',DR) 
300  RF:nt)<'i,3iO>     liDR^R.O,  <TRLKmD/.T>.DRTn<IiJD,  J>,  J=I  ,m)R) 

3i0  F0K(1Rr<R2,62RCO 

DO    -^00    LP3^1,36 


DRP0<i>=DRTR<i>LP3> 
RPCi<I>  =  G.O 
KP  - 1 
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DO  500  i.r4==i:po*5 

KI-"»KP+1 
jm)^LF'4+l 

I)f;rO<KP>  =  DRTn<lHD,LP3> 
&eO  RPCi<KP)=LP4 

CJ'S.G 
C.V«  1.0 

Fl-6.  G*<DRPCi<3)-t)RP0<2)-<t)RP0<2>-DP.PCKi  >>''3.  6> 

F 2-^ 6 .  &^' <  DRPO  <  'i  > -ORPO  <  o  )  -DflPO  <  3  )  +1>RP0  < 2  )  > 

Pj^=C'.0 

P2"<F1*C4-F2)/'CC1*C4-C3> 

P3"Fi-CUP2 

P't^C.O 

DC    cUtC    LP5-1.-MP0 

nr,:=^LP5+&  - 

P(i^;=^rL0RT<l.P5> 

DUJU"'-P2/IS.0*P03>S'i:3+<DRP0<:2)^3.6-P2/2.e)*P0S 

I     +<DRPoa  >/'3.  e>*<:3.o-pos) 
ir<Duno>eict.6t'-6,6DC 

6K<  DLilfO=0.  D 

GOG  DR7R<IMD.LP3)  =  DLIM0 

'IGG  CCUITIMUE 

DC    ?00    LP€.=  1.|{P0 

l..'RJITE<3>7C2>    G 
?02  FDRKRT':'!',  iCX''Q='F5.  !///'/'> 

DC    ?C0    LP7=2,2 

LRC==5*<LP7-1> 

DO    7'50    LP8=1  ,3G 

)Hi>2^LP6+l-{LRC 
DUin-4)RTR<ilU)2/LPS> 

Dcno-iGC.  O'^Durn 

l'.Ui;G^-lFlX<DUnO> 

I  r<i5UDG-5CO?GS.  705,706 
7GC  HUl:G-9e 

701.  DUI'iG=^G.5H.Fl.0RT<LP8) 

1F<LP8~1  >7'^5/?5r.,756 

7f.r.  e.ui^  soM-sDuin 

GO   TO  ?r.G 
?f.f.  J  F  <  i.  Pi!- 5G  >  7T.7 ,  7f.5  >  75D 

7t.7  £;UI^=SLI|U?.G^[>UK1 


11^ 


?10  FCiKi;Rr<F?.2/E20.  ?,5.X>  lOCtftl  > 

l'R]TL<3.?2tO    F'ROB 
?20  FCiRr5RT</'/'-'lCiX''PRC!l3fin:lLITV=''F7.4/'/'/'/'> 


£: 

CLH 

CLL 

6 

TRD 

>LF^6 

£ 

TRf> 

<6e 

f. 

F.'TL 

FnLiftTL 

c; 

Tftf: 

<56 

s 

DCfl 

\[)R 

S 

CLU 

S 

TftD 

\LP1 

c. 

TRfJ 

<<;.C' 

S 

RTI. 

RTLiFaL 

S 

t>Cfi 

\DR« 

fi=^rLCiftT<LP6> 
Cfll-L    OOPEIK'SVS'.DR) 

IMHTFC''.  ,?91>    K[)R,R.Q,  <TRUC1  ,  J  )  /  DRTR<mD2.  J  >  ,  J=--I  /|{l)fi> 
?9i  FCir,l1R7(R2^e.2P.6> 

CRLL    OCLOSE 


?C'.t< 

CCiMTIKUE 

2(JCV 

.  ccirniMuc 

EHIi 

115 


Pf^OGRftM    IHPOL 


PI  nrN's I  on  Tnuc 3 1  > ,  Dfiin  <3 1  > .  mug <  i  co > ..  ptr c  i  se > .  PDfK  ksc? ) ,  cci ,  -i  > 

t>ll'.ENSlCiU    F<'i;>I)f;P0<6>,P<:6) 

[>o  r.o  1  =  1 ,  IOC 
UD  nuo<i>=-i3?& 

C:<:f/I>=1.0 
C:<2i3)=I.O 

c:<?.  n=e.o 

C:<3.2>=^I.C  -  ^ 

C:«:?.3)  =  4.G 

C<'5,3>  =  1.0 

Cni.L    l5ftTllKC,!1.D> 

DO    108    LP I =i, 5 

DO    100    LP2-3.9 

S  CLP.    CLL 

S  TRD    \LP1 

?.  TP.D    C60 

£  RTL;RTL;KTL 

S  TRP-    <56 

S  DCf!    NDH 

S  CLL 

S  Tnr-    \LP2 

S  TRP    (60 

1^  RTL;RTL;RTL 

S  DCH    \bm\ 

19?  CRLL    10PE:K<''S:VS',DR) 

RrjP.DC'l,  I&CO    W[iFl,n. a.  <TRLKI  >,DflTR<l  ),  1=2,31  > 
if^-e  F0!d1RT<R2,t-.2R6> 

TRli(l>  =  e..O 

DRTR<l>=e.0 

KHKP^O 
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DO    ?t^0    LP3-=I,26 

im5=^LP3 

[>0    300    LP4=I/6 

Df-:f  0<LP4>  =  [>RTn<IH0> 

DO    326    I  ==1.4 
32C'.  F<)  )--6.  e.*<[)nPCKl+2)-2.  =5t>RP0<  I  +  I  >  +  DRPO<:i  >> 

DO    350    }~1 ,A  '  ' 

PC)N[)>  =  »3.0 
DO  350  J==I,'1 
356      P<3!il))  =  P<IHD>  +  C<l.  J>*F<J? 

p<5>=^e..o 

P<(.>=6.0 

1F<LP3-I>410..  4  1G.40I3 
4(tC<  ]F<l.P3-26)42tS/43e*430  ' 

4  JO  1.00  =  0 

GO    TO    450 
430  L00=--2  • 

CO    TO    450 

420  JJ{fii"!:EEP  I 
DO    421    LF'5=l/5 

DTfill=0.2*Fl.0RT<LP5> 
PTfi<)IU)l>  =  0.  i>^FL0RT'Cl(HDl> 

PI>R<IKC>1  >  =  P<3>*<1  .  0-DTfiiJ?*-i3/'6.  0  +  PC4  )*DTRU**3-^6 .  0 
1       +XDRPCK4  )-P<4  >/-£..  6)*t)TRU  fr  CDRP0<3)-P<3)/6.  0)*<  i.  0-1>TRIO 
H  <PDRC]l<t)i>>425/425,42i 
425  P<3>=0.9«P<3> 

p<4>=e.9*P<3) 
GO    TO    426 

421  coin  1  HUE 

iu:i;p=iNDi 

GO    TO    200 

450  1»!)1=KEEP 

DO    451    LPe--l  ,3 
DO    451    LP7=i>5 
IWDl'-IUDl  +  l 
DTnU=--0.  2*FL0RTCLP7) 
PTR<1HDO=:0.  1>:^FL0RT<1HD1> 
K--^I.PC-.+l.OC 

f'DRClMDn  =  PCK>*Cl  .  0-DTnU>*ri3.''C. .  0  +  P<K+l  >*C>TR0**3/6,0 
1       -4  <OliPD<K+I  >-P<K+l  >yS.  0)'«=DTnLl+<DRPCKiO-P<l<)/'6.  e>--«<i.  8-DTRU> 

451  coirniiuE 

KEKP-lMDi 
GO    TO    200 

200  coininuE 
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sun  :=  0.0 

DC    C.GO    1=1,150 

PKOE:=SIJn*0.  1 
PO    C2ti    1=^1,  I 'JO 
d.20  PPP<]  >^p[>fi<I>^PROB 

KJ)f;=^i50 

CfiLl.    00PEK'<'SVS',[)R> 

IM::1TF<4,630)    MDf-|,R,0,  <PTn<I>,PC»R<I>.- I  =  1/HDR> 
f.?0  F0rinRT<R;?,362R6) 

CRLL    OCLOSE 

100      COHTinUE 


il8 


F'RCigrrh  pco»3 


D  J  1U:K'S  1  on    f "DET  < 20)  ,  T L  < 26 >  .  TM < 2CO  ,  T  1  »U < 20 >  ,  SMRL  < 20 )  ,  KRT  < 2Ct ) 
t>ini:i{S10l<    7RIJ<I50)/DRC<<U.e)/DHl<lliO),MUG<160) 

DO    DC    1=1,300 
•oP  KUG<1>  =  -13?6 

CRLL    lOPtKC'SS'S'/ 'ZXK'O') 

Rfif^X'i-c.o)  iu>n,  <:Ttti.i<i>,t)RO<i ),  i-i  ,iipn) 

f.O  FDRnRT':ft2.302R6) 

DO    IDG    LPI=I,5 

t>C    100    LF'2^3/9 

S  CLR    CLL 

S  7Rr>    XLPI 

f:  7RP    (60 

S  RTI.  iRTLjRTl. 

t;  TRr:-    <:t".6 

S  DCR    NDR 

i>  CLL 

J::  TRr>    \LP2 

S  TRr>    <60 

S  RTLiRTLiRTL 

S  PCR    NDR4t 

J9?  CfiLL    10PE:H<''SVS</[>ft> 

Rf{RD<4,26G>    «f)fl. R^a,  CtRlKI),[)Rl<l>,  1  =  1, NDR> 
2^0  F0RI?R1<R2,3O2R6) 

GO    TO    220 
l'RiTE<:3,210>    R,Q 
2J0  FOM'.RT':'-!  'lOK..  'R^'FS.  6,  IOX>  'ft='F4.  !>'//'> 

220  J=0 

DO    400    LP3=1,2G 

J==^  +  I 

PDHKJJ^l  .  0-0.  605*FL0R7<LP3> 

7P«O0.0:<PDE7<0) 

HOLD^DRKO 

LOL^-1 

L0I-!=1 

1  WD  1=^1 

£F.RL<J>  =  1.0 
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DO    360    LP4=I,149 

K0L0=M0LD-0.r.*<[)fiI<L0L>-fDftKLP4>) 

L0L=LP4 

sun  1= HOLD 

3  jr.  H<SIU'.I-TP>32P>330*35O 

326  J»!ii-llU)l  +  I 

1  P  <  1  niy  i  - 1  50  >  322.  322 .  360 
322  HOLD=SUiU 

S:UIU-SLIin+C.5*<DfiI<l.CH>  +  r>«i<:iWD!>> 

LOH^IHDl 

GO    TO    3iS 


336 

sun 6= 0.0 

I>0    33^.    1^L0L..10H 

1  F  (  i  - LOL  > 33t-:. ..  336 ,  337 

33? 

1 F  <  I  -L0n>33r:! ,  336  ,  336 

336 

sun 0=SUI 50+0.  ti*DftO<  I  > 

GO    TO    33S 

33f^ 

SUn8=-SUI'.04DR0<I> 

33'o 

coininui: 

PR0E:=SUH0*e.  1 

Tlh'T^-O.O 

GO    TO    366 

sf.e  t>ur<2=sut'.  i  -hold 

DUin==TP-HOLD 

Tiia=t)Uin/'pui52 


sun 6= 0.0 

LOi-:  =  l.OH-l 

DO    3t^P..    1==  LOL..  I- OH 

]  F  <  I  -LOL  >  356 , 356 , 357 

3'o? 

I F  f 1 - L  OhO  35S , 356 , 356 

31-6 

SUI';C:^SLii;0  +  0.  5*DR0<  I  > 

GO    TO    355 

35f! 

SUi:6^-SUI*.0+DR0<I> 

3'3r. 

coiainuE 

HOLDO-^SUItO 

SUno-SU!  16+ 0 .  5*  < DRO  <  LOH  >  +  D R0  <  I  l{D  I  >  ) 

PROB-::  <  MOLDO+  T  I  WT >;••<:  SUiiO-HOL  DG  )  >  ^-Q .  i 

IMPi^^LOn 
360  ]F<Pr<0B-SMRL<J>>376/ 360,300 

370  Snf-;L<J>  =  PR0B 

tl<:j>=tru(lol> 

TH<J)=TRUCLOH>+TIin*0. 1 
T  J  IR« <  J  >  ^T  H  <  J  )  -  T RU  < L OL  > 

360      COHTIHUE 
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RftT<J)  =  PDETCJ>/'SMRL<J) 

60    TCI    A&d 
596  I'Kl  TE  <3,  392  >    POET  <  J  >  ,  TL  <  J  >  ,  TH<  J  )  ,  T  1  »U<  J  >  ,  SHRL  <..l  >  ,  RfiT  <  J  > 

392  F0rvHftT<'1F€:.3,2E15.  ?> 

100 


C.Cl|^• 

TIIRIE 

K'DR 

-J 

CLM 

CLL 

TRH- 

\LF'I 

TRP 

<6G 

RTL 

!fnL;RTL 

TRP 

<.46 

[iCn 

\C)fi 

CLl. 

TRP 

\LP2 

TRP 

<i.{i 

RTL. 

RTLiKTE 

DCfi 

xDfii? 

CP.LL    OOPEJU'-SVS'.Dft) 
t'RlTE<4,mO)    NDR, 
1       XPPETC  I  ^  ,  TL  CI  >  .  TH<  I  >:  7  I  HM<  1  >  ,  SMRL  C  1  >  ,  RRT<  I  >  ,  1  r:  i  ,  |^|>ft) 
1J0  F0r;;i;RT<n2,260B6) 

CRLL    OCLCiSE 

JCtP  CCUniNUE 

EMU 
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PROGRflM   ftCO» 


t)3!1EHSlCif<    PD<2O),TLC2e),TH<2e),TI<?6>,PF<20),Rfl<2CO,F<2n 


F<}  >  = 
F<2)  = 
F<?)  = 
F<'1>  = 
F<fO  = 
F<C.>  = 
F<7)  = 
F<i!>  = 

F<:^;>= 

F<UO 

F<n> 

F<J2> 
F<J3) 
F<i4> 
F<lli> 
F<  J  C > 
F<J?> 
F<HO 
FCJTO 
F<2C) 
F  <2 1  ) 


l.Cf 
1.0 

2.e 

24.8 
126.  C 
720.  0 
50-4  Ci.O 
'^•3326.  0 
=  3.fc.2SS*l 
=  3.  t.2fc:S*l 
=  3.99168='; 
=4. ?9f la. 
---6.22?>iU« 
=  S.?l?e3^: 
=  i.  36?e.?* 
=  2.  6922$:=^ 
-3.5568?* 
=^^6.  -18 23? I- 
=  I.2i6'^5* 
=  2.  4329e.'N 


10. D**? 

io.e**io 
io.e*=;:i2 

I6.0+i^I3 
18.G-i-^14 
10.0**15 
10,0*'^^1? 
10. 0**18 


[)0  ICS  LF"i==4.5 

DO  ICO  LP2=4,4 

CLfl    CLL 
TR't?    vLPi 
Tf5r>    <68 
RTLiRTLiRTL 

[>Cr^    \Dft 

CLL 

TRr>    \LP2 

ffm  <60 

PiTLiRTL/RJL,  , 
t>C:fi    Nr.)f54i 


no 


CfiLL    lOPF.KK'SVS'.DR) 

RKn[)<4,110)    KDR/ 

CPIXI  >,7L<l>,fH<:i),TI<I>,PF<l>,Rft<I  >,  l=^l>NDft> 

F01:nRT<:fi2,200B6) 


l&O 


0"(<.  I'JrL0RT<LP2> 
lMi]TL<3*  120>  LPl.G! 
F OM'.RT  <  '  i  '  1  OX  "  R=  "  1  2 . 


10X'0=^'F4.  I///') 
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[)0    iCC    LF'3=4,8*2 
H«2*LP3 

DO    100    LFM«IO, 10 
fiL"6.05+FL0ftT<LP4> 

DO    J  00    i-F-5=Lr'3iH 
t»KlTr<3.  130>    LP3,ftL.LP5 
I  30  FOUUfn  <  i  X ,  I  1  0  >  F 1  0 .  2 .  I  1  0  ) 

DO  300  LPO=:i  ,18 

suni^o.e 

UK-n-LF'^^l 

DC  260  I.P6=I,I:H 

C"rCI-lH>^F<LPt.>/'F<N-LP642> 

DL.'I'0=PF<LP&) 

ni!i;2=--LPt.-i 

200      SUin--=SUin  +r.*DUt16H'i=MUDI  *<I  .  6-DUM0)**fRID2 

60  TO  226 

IFCSUIU-e.  00^".>  100,  21  0,210 
210  JF<i;.UIVl--nL>22e,220.-300 

220  f;UIT2^0.0  , 

DO    250    LP7--=i  ,KM 
C^F<HH:i/'F<LP7>/'F(H-LP7  +  2> 
DUnO=PtKLPS> 
Kl!|;I=^H-LP7+i 
r:UL:2-LP7-I 
250  SUn2==Slin2+C*DUI1O**HUDl';:<i.0-DllHO>*«HU&2 

Dl!!<0:=SUH2>'SLli'il 

1-'F:1TF<3,270>  LPS,PD<LP8>,PF<:LPS>,SUE'.l,SIJt12,DU!tO 
270      FOkHRKIOa,  13,5EIC..  ?) 

300      COKTIHUE 
100      CONTINUE 

END 
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PROGRl'lt!   F-LIM 


PI:  3.  I'nr.92? 

ni=  D.  3 1938 IS 
f{-?=  -0.  3f.t.563S 

ft.1^-1  .821256 
fi;:^=  n  336274 
LfiV=^10 

DO  100  l-F-i  =  l  ,3 
R"rLORr<LPi> 


ijct 


t)D    IOC    LP2^4.'? 

P.-e.  l^=FLafiT<L.P2>  , 

1'}-.:itl:<3.  110)  LPi,Q 

FDRHrtK-"!  ''I0X''fl='I2.  IOX'C4='^F4.  l/'x'/') 


1  :>fr. 


r>Ci    100    LP3=4,8.2 

T  =-2 .  0*P I  *FLOiir  <LP3> -^Q 

l'f-::jTF<3,  120)    LP3 

F OR'-JRT  <  X5X ' »i=:^  ^  1 2 > 


sun  2-.^  0.0 

I>0    260    LP4  =  l.LRy 


220 


230 


200 


c-.ur. 

K-K 

H  < 
DU/' 
S :;  f ; 

mm 

GO 

sun 


P4-1 

L ORT  <  HUD  )  *P i  ^'F L ORT  <  L  R V ) 

0.0 


t-:  +  FL.0RT<iO*PI 
■T>23S,230.260 

<DUUO) 
SUni+S*'?-2 

220 
SUIt2+SUiU 
5*SLIlt2*R**2-^FL0RT<LRU> 


DO    100    LP5=I,26 
KU!:=^l.P5-6 

[>unor:  j .  o=<-FLORT  <:mud> 

UH  I-  -  DUnC*  RBS  <  D  UMO  )  /-P I 


12i|- 


>{C=-<tlf»UEL>/<2.G*EH)**6.5 

s>-- 1 .  K/- <  1 . 0+  Ri?.*t)un  1  > 

PC*.:- ?xr  < ft  1  -:  V■^  ri2 V  Vw  ■<:2+ R3 *V*'«  3  +  fi4 1=  V* *4  +  fi5* V+*5 ) 
ir<XrO2t.t:,2D0/251 
25C'.  F'O"  I  .  C-F'6 

251  >:i==C--EM+UtL>/'<2.  D+EfO+^e.S 

l>Ui1C;--Xi»;*2/'2.e 
2- 1:  J^f  < Dunct )  /-  <  2 .  6>t  P I  ) **i3 .  5 

F*  1  --  r-^  >»:  <  fi  J  ^  \'  f  R  2  «  V  *  i-  2 + ft  ?  *  V  *  :i-  3 + f  1 4  *  S'  *  ^  4  +  ft  5  *  V  *  *  5  ) 

JFCXi>2&ft, 2601,261  _     „_  _  

2C0  P  1^=1.0 -Pi 

2C.1  RRT==Pi/'Pe 

l?R  1  T r;  < 3  /  360  )    L  P5  >  UEL  ,  XO ,  X 1 ,  PS  i  P 1  .  RfiT 

SC'.C.  F  Oft  (IftI  <  1  OX  /  13 ,  6E  1  5 ,  ?  ) 


100 


COf'lTIHUE 
El^!> 
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PROGRRfl   F'SEG 


PJ"3. 1115927 
RC<=-0.23I£.4I9 

R3'^I.?e.M7S 
R'J- -1.821256 
RC" 1 , 33&274 

DO  100  LF-S=2,  160,2 
C.'.l<==0..  01»;FL0Rf  (LPS) 

00    100    LPI==I,2 
R-fLORTCLPl) 


no 


DO    100    LP7- 1,11 

i;ut;^LP?-i 

RI,'"R*FLOfiT<IJ!JJ>) 


DO    100    LP2=4,4 
G==0.  l>;FLaRT<LP2)  , 

l!ri:lTE<3,110>    R,D,RH,OH 

FORHRT  <-''/','//-' yxi  OX''fl=''Fr..  2,  IGX'Q-'FS.  2j  lOX'RM-TS.  2 
1        ,  10X't!!{-^'F5.  2/'''-''>  I 


1  20 


DO    iOO    LP3-4,t.,2 
T«=2.  0+Pl';FL0aT<LP3>/'e 
t'RnE<3,  126>    LP3 
F0Ri1RT<.''5X'N-=''12> 


?Ui^2'^0.O 

TSUKH=0.e 

DD    200    LP«1=1  ..LRU 

I1UI>  =  LP-1-I 

T  H\:  ^  F L  ORT  <  riUD  >  *P  I  /F L  ORT  <  L  flU  ) 

$:IUU'^0.0 

SUKK^O.O 

K=-l 
220  K--K+I 

Ti:==TMf-:  +  Fl.ORT<K)M-PI 

IF  <  TIC-  T  ?230  ,  230,  21  0 
230  DUnO^C'-i'TK 

f.UItl-SUI'.l-f  s>?=;-2 

DL»f:G:=aN'>JrK 

&N^=^SI{KI)UIS0> 

GO    TO    220 

210  Eutn2-^sun24suni 

200  .  Ti:UnK'=TSLI!1!M  SUI'lH 

F(<"0,  r.*SlJI(2s.n**2-^FL0nT<:LRU) 
t:H2-TSLiHH'frt>:.n.W/Fl.0fiT<l.nV>--t:H 
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PD  100  LP5=1,I5 

I)Ur.O-l.e:>iFLDRT<Mll») 
UH  I. '  DUKO^  RBS  < DUHO  )  ^P I 


260 
2fci 
300 
100 


I)lM'0--Xe.H*2/'2.  0 

2~  i:  y.T-  <  DUI 50  )  ^-  <  2 .  e*p  I  >  *« 0 .  5 

I)Ui:l=:FiBS:<XO) 

V"  J  .C.''<:I.C-^RO^DIJ^I) 

lF<X0>2i:.U,2r.0,251 

P0:-1.C-P0 

[>Ur.O^--X  1^*2/2.  0 

2~iiy.f'<.t>imoy^<2.  6*pi  >**0.5 

DUf-l  =  nE:S':>;i> 

V"J  .0.''Cl.G+r{C«DUMI> 

PI  ^-ifi^R  i  ^:V+ft2>:  V++2+ft3*V**3+R4*V^*4  +  ft5HV**5> 

H  <Xi  >2i:.0,2fc.8,261 

Pi=-1.0-PI 


RRT^PI/PO 

rRitf:<3.3SC)  XP5.VEL,X0,KI.PG,PJ,RRT 

FORMnf  <  I  OX.  1 3,  t.t:  lli .  ?> 

cctsrriMiiE 
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PRD&fiftM   PCOtH 


COniJOM    TftU/DRC.DRl/MUG 

I) I !'.£:!{•:": lOM    r-pt:T<:2e),TL<20>.TH<26),TllKK2e),SMP.L<26),RRT<2O) 

DO  r.o  1  =  1 ,  100 

GO  I'.UC:<1>  =  -13?6 

i?irjTi;<i  ,57> 

5?  FOFdlRKlX'TVPE    IH    THE    REFEREf-JCE    DRTn.') 

R[rHD<:i.S9>    DRS 
59  f  Dl"vMnT(ft»S> 

CRl.l.     lOPENK'SVS'^DRS) 

RERIX-1.60)    NDR.RS.OS,  <TniKI  >/DRl  <3  >,  1^1  ,HDR) 
(-.0  F0KI«RT<R2,302ft6) 

[>0    iOO    LPI=2,2  »  ' 

DC!    100    LP2=3/9 

S  CLR    CLL 

S  TRP    \LP1 

f:  TRD    <€.e 

S  RTl.  ;F:TL-RTL 

S  TRP    <5& 

S      .  [)CR    NDR 

S  CLL 

f:  TRr>    NLP2 

t:  TRr^    <6e 

£  RTLiRTLiFaL 

S  t>CB    \[>R1f 

199  CRll.    lOPEfJC'SVS'^.DR) 

FvF:ri[K4.2C'iO)    Nl)R,RH,QH/ <7RU<I>,[>R0<I>>  I  =  1..HDR) 

200  FDKMRT<R2.3C-i2R6) 

l=!FaTE<?,210)    RS,OS.ftN,QH 
210  FOi^l'.RK'i'lOX,  'RS^^'FS.O.  16X,  'QS=''F'i.  I  ,  1  0X'fiH='F3.  0,  1  OX 

1       '(.'!{= 'F4.  !/'/'>'> 

220  J--=e 

t>0    400    LP3=I,20 

J--J  +  i 

Pl>f:T<J>  =  i  .  0-0.  05*FL0RT<LP3> 

TP=10.e';PPET<J> 

HOL[):^^[)Ri<l> 

LOL=i 

L0!>=1 

INlU'^l 

Sr.FL  <.!>:=  1.0 
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I>0    300    LF'4  =  1..  149 

HOLD=HOLD-0.  r.*<Dni  <L0L)  +  t>ftI<LP4>> 
L0I-  =  LF'4 

3  HV  1  F  \  SU! !  i  -  T  P  >  ?26 ,  338  /  350 

326  1MJM=--I»<C>1  +  } 

1FC1MD1-U.CO322.322.3G0 
322  HDlD=:SLIi11 

SllJ'l!:SUi;i+G.  lf.*<Dfll<LOH>  +  Dfil  <IHC)1  >> 
GO    TO    :sUi 


330 

SU!'0=0.0 

DO  335  l^LOL.LOH 

I  f  i:  1  - L  OL  >  336 ..  336 ,  33? 

33? 

1  F  <  1  -•  L  OH  >  330 ,  336 ,  336 

336 

SUfTC5=--&UitCitO.  t'.*Df-iO<  I  > 

60  TO  335 

330 

OllJ'0-SL!nO+[)ftO<  I  > 

33f. 

coj-niijut; 

FRrjv.=i.mi(i*Z<.i 

Tiin^-o.ci 

GO  TO  369 

315.0  l>UK2=^SI.I(n-H0LD 

DUI' J  =TP~ HOLD 
T I  HT--^[5U|{  1/00152 

£UriGr:0.0 

LOIi^-LOH-l 

[)0    355    1=:LOL,LOH 

1  F  < 1 -IDL  >  356 , 356 , 357 
35?  1F<I--L0H>  358. .356,356 
356      6lU'.0"SUIiO+i;v.  5*DnO<I> 

GO    TO    355 
350  SOI'iO~OlJ;iO-i-DFsO<I  > 

355  COJaiilUE 

HOLDO^-OUno 

suiio===suitcs-to.  5*<:t>fiy<:LOH>+[jno<  I WDi  >  > 

PKOO^  <  HOL  [J0+  T  1  HI  ^-^  <  SUMO  -  HOL  DO  )  >  ■♦=0 .  I 

in(>i=:LOH 

360      IF<PF;0E!-SHRL<J>)3?e, 300,300 
3?0      Sr'!r;L<J>:=PF'OB 

Ti.  <:J):^TftL!<LOL> 

T H  •:  »i  >  -^T  nu  <  1.  OH  >  +  T I  in H  0 .  i 

TlKU<sT)~TH<.O-TflLKL0L> 


300 


COiniSiUE 
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KflT<J)=POET<J>/'SMftL<J> 
390  lM.:irE<3,392>    PDEK  J>  ,  TL<  J>  ,  TH<  J>  ,  T  iHM<  J)  ,  SMRL  <  J  )  ,RRT<  J> 

392  FCiRimT(:'iF8.3/2EI5.?) 

40Ct 


coin  I  HUE 

KI^R 

^J 

cin 

CLL 

inrr- 

NLPl 

TRP 

<€.6 

R7L 

RTliftTl. 

THP 

<46 

t)Cfi 

N[)ft 

CL.L 

TRD 

NLP2 

TP.rs 

<6e. 

RTL; 

RTL;RTL 

DCfi 

\i>M 

CfilL    OCiPEK'<'SVS',Dfi> 
I'RiTEC'^^'iUt)    \U>B, 
i       <f'!>ETa>,TL<J  ).TH<I>,TllRKI)-SHfJL<I>.RfiT<I>,  !=-!  ,W-I>ft) 
'IJO  F0r::MnT<fi2,2e.0?H6) 

Cfii.L    CiGLOSE 

3  GO  coin  1  HUE 

EH15 


130 


PRCiGRrjH   F'K'OI 


PI-3.J '11592? 
P.tu.0.2?l.;.'1I9 

fi2-  --0.3^-t.::.638 

n5--  I  .  33Ci274 

DO    IGO    LPi=^I,3 
R-J-LORTCLPl) 

t)C    iCQ    LP7=^2,IS,2 
CH0I=FL.0nT<LP7> 


JJP 


[)D    J  00    LP2=3.9.2 
a=-C<.  l*rL0r-lT<LP2> 

t'RrrE<3,nco  fi.o.cNoi 


3  2e 


t>0    100    LP3  =  -4,8.2 

T'-:^.0';pi  ;FLorn  ■cLP3>/'a 

l'RITF<3>  120>    LP3 


sur,2--ci.o 

DO    260    LP4=I.LRU 

irLU>=Lp4-i 

1  MM  =  FL.ORT  <  KUD  )  *P  I  -'PL  ORT  < L RU  ) 
S;UIU==0.0 

220  \>\:n 

..   TK=^TH(r  +  FLORT<:K>*PI 
1F<T!:-T>230>238,2e0 
23^.      Di:i<0=C'-fTK 

s=-i;j)ur)UHO) 

«:UnS:-S:!JMl4S*?^2 
60    TO    220 
200  Slll12^^SU^2•^SLIIti 

EI^^O.  f.*SUn2>:=fi**2/'FL0RT<LRU> 


DD    100    LP5=1 . 10 

DUHO^- 1  .  OH-FLORT  CMUD) 
V{:L-=L)Ui!0**2^PI 
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if5C. 
251 


C)ur:C:r.-xoi-v2/'2.e 
z=^i;Kr<:[)ijnco/<2. 0+pi>**0.5 

I)in',ir.-RE:S<XO> 

rev-  v^.  <  fi  1  *  V  ^  ft:?*  V+  +  2+ ft3* V* >»  3+ R  4  -s- V>?  * 4  ^  ft5* V**5  ) 

IF<XO!>2?:*a,250,251 

PC'.-  I.Ci-P9 

>:J--<-E!!+UEL)/'<2.Di:EtO*«0.5 
[)lli'.C<^-Xl**2/2.  6 

r>uni:^fiF:s;<>;3  >  

s'"  I .  (?/<  1 .  o+fte^cJLirii  > 

F•3'-■?-^<RHY•^f{2>:  V*-i:2+fi3*V*-'<3+R4'«:V**'l+ft5*V*-'<=5> 

iF<Xi>2';.G..2t.e..  261 

P1:U.G-P1 

RRT^PJ^pe 

I'f? }  T E  < 3  i  380  >    LP5  A'EL  ,  NO ..  X  i .  PO ,  P 1  ,  RRT 

FOFd1R1<iOX..  i3,6E15.  ?) 

c-oininuE 

EN[5 
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